jlperla / continuous_time_methods

MIT License
23 stars 15 forks source link

Solve by multiplying by $\Delta$ rather than diving $A$ by $\Delta$ #24

Closed jlperla closed 7 years ago

jlperla commented 7 years ago

You will see that I modified /matlab/lib/simple_HJBE_discretized_univariate.m (in 76e937325bbd50720128a46ea2f20931e7ac9cbb ) to NOT divide A by $\Delta$. Instead, it solves the system with the $\rho $ and the utility mulitipled by $\Delta$. This is as written in operator_discretization_finite_differences.pdf equation (39). The purpose of multiplying everything by $\Delta$ is primarily numerical stability for small $\Delta$.

I didn't make the corresponding changes to simple_optimal_stopping_diffusion.m or simple_joint_HJBE_stationary_distribution_univariate.m, but we need to do a similar change to avoid dividing by $\Delta$ in simple_joint_HJBE_stationary_distribution_univariate.m and simple_joint_HJBE_stationary_distribution_univariate.m and make sure that the test cases pass (storing off new versions of the files as required).

stevenzhangdx commented 7 years ago

I have modified the other two files and all tests passed after modification.

jlperla commented 7 years ago

I am reopening this because it isn't what I was looking for. You will notice that it NEVER divides by $\Delta$ in the implementation of https://github.com/econtoolkit/continuous_time_methods/commit/76e937325bbd50720128a46ea2f20931e7ac9cbb

The code v = (Delta * rho * speye(I) - A)\(Delta * u) isn't really dividing by $\Delta$, but rather solving the linear system: $(\Delta \rho I- A)v = \Delta u$ which is a subtle but important difference.

In simple_optimal_stopping_diffusion implementation in line 71 you are dividing by Delta. What we want to do is some transformation of $q \to \Delta q$ (which may change something else in the math when you go through the transformation of the problem.

For example, go to the optimal_stopping.pdf notes and make sure we can implement the problem exactly. Is it possible that (5) is incorrect? Perhaps it instead be $ 0 = \min{\Delta\rho v - \Delta u - A v, \Delta v - \Delta S} $

Similarly, for simple_joint_HJBE_stationary_distribution_univariate.m it is dividing by Delta in line 44, but I want to solve the formulas exactly as specified in equation (41) and (42) (or show that they are incorrect if!)

One thing that helps is putting the exact line number from the pdf for the code. If you need to divide by Delta or something like that, then maybe the PDF needs to change.

stevenzhangdx commented 7 years ago

I have removed all divisions of $\Delta$ and followed equations in the notes.

jlperla commented 7 years ago

Looks great. I noticed a few more $\Delta$ divisions (i.e. in stationary_distribution_discretized_univariate.m) but played around with it and also fixed the document to reflect that the LLS needs to change slightly. You can get the latest to see what I did.

I also noticed that the testsuite was failing on `'simple_optimal_stopping_diffusion_test/baseline_HACT_test' which looks like it may have been due to an overwritten testfile. I updated this file and added it to the test. The precision of the comparison was decreased since this one came from the original dividing by Delta code, which was numerically less stable (since it was /1000^2 or something like that). This was the reason to get rid of the /Delta originally, so not surprised they don't quite match.

stevenzhangdx commented 7 years ago

I read what you did to LLS part. But I think there is no need to divide $\Delta$ in (31), because in my understanding multiplying every element in a matrix by the same constant only rescale eigenvalues but eigenvectors. It is not a big deal since the result does not really change.

jlperla commented 7 years ago

Thanks for checking on it. Yes, you are correct that we don't need to divide by $\Delta$ in (31) to find the $f$ in practice (even if it is the 'correct' formula). You will see at the top of page 6 that it says we multiple by $\Delta$ to just find the eigenvector for the $A$. That way we never need to divide by an extra $\Delta$, which is the key goal. I think that is how the code is right now.

Now, with the LLS part it is a little trickier since it isn't solving as an eigenvalue problem. In that case, since the normalization is built into the system of equations we need to be more careful with the $\Delta$. You can check to see if (33) to (36) seems correct.