Code and test cases for generating a dynamic operator

[jlperla]

After #31 we can generate the $A$ matrix and verify it.

I added in two files \lib\discretize_time_varying_univariate_diffusion.m has a skeleton of the function to generate the A, just like in the non-time-varying case. and the "matlab\tests\discretize_nonuniform_univariate_diffusion_test.m" provides a skeleton for it (including the way to stack the $\mu(t,x)$ and $\sigma^2(t,x)$ correctly).

Once you generate the code for the $A$, we should add in all sorts of tests to make sure t is correct.

[jlperla]

For this, lets start with a few very simple dynamic examples

• [x] First, use a time-invariant $u(x) = e^x$ along with GBM $\mu \equiv \mu x$, and $\sigma^2(x) = \sigma^2 x^2$ to make sure that the method correctly finds a time-invariant solution to $v(t,x)$, and that this matches the stationary solution we solved separately
• [x] Next, make the utility slightly time varying with some parameter. Maybe $u(t,x) = e^x F(t)$ where $F(T) = 1$ and F(0)$ is parameter by something to play around with
• [x] Then leave $u(x) = e^x$ time invariant, and have the $\mu(t,x) = \mu x F(t)$, again where $F(T) = 1$ and you can mess around with it.
• [x] Make sure it works (or at least that the code is checking for regressions) with non-uniform grid points in both directions.

For the $F(t)$ function, consider starting with a linear function passing through two points $(t,A)$ and $(T,1)$ I think the formula for this something like $$ F(t;a) = \frac{T-a}{T}t + a $$

For these, we should play around with time step sizes (and whether it is uniform or not) to check if the method is stable. If there are major stability issues, it suggests we really need to use the implicit setup.

[jlperla]

I came across this: http://math.mit.edu/~stevenj/18.303/lecture-10.pdf

It is for 2 dimensions, so not directly applicable, but has a good explanation for playing around with Kronecker products for generating the finite differences cleanly.

[sevhou]

I've revised the discretize time varying diffusion file in \lib path.

The code is a non-uniformed one and can be easily applied to uniform case I think. The code is for explicit time step only, where $h{+,n}=t{n+1}-t_{n}$.

For the tests I haven't finish the coding. But I manually checked a small $I$ case which seems work fine. Will continue working on point 1 and 2.

[sevhou]

Guys I've finished point 1 in the test, and a basic test for time-varying method.

The basic test tests a time varying $u(t,x)$, $\mu_{t,x}$ and $\sigma^{2}(t,x)$ case and use the time-varying method to solve it. The code works fine and the permutations work the way they should be.

The \emph{non-time-varying-test} tests $u(x)=e^{x}$, $\mu(x)=\mu x$, $\sigma^{2}(x)=\sigma^{2}x^{2}$. The code solve the same system using time varying solver and non-time-varying solver. It does the following comparisons:

1. Compare the sub matrix $A_{1}$ in time varying result with $A$ in non-time-varying result. The difference lies smaller than $1e-5$.

2. Compare the sub matrix $A{1}$ and $A{n}$ ($n$ randomly drawn) and the differnce lies smaller than $1e-9$.

3. Compare the value function from time varying result $v$ with that from non-time varying result $v_{non}$, difference lies smaller than $1e-5$.

4. Compare the value function when $t=h{1}$ and $t=h{n}$, because the $v$ we get from time-varying code is $R^{NI}$, the difference lies smaller than $1e-9$.

Will continue to do the other tests in the rest of week.

[jlperla]

See my change to the code to make the A sparse. I had trouble running the unit test, but I don't think that had anything to do with the code I changed.

[sevhou]

Added time-varying-utility-test, which does point 2 basically: $u(t,x)=e^{x}F(t)$, is time varying, $\mu(t,x)=\mu x$ and $\sigma^{2}(t,x)=\sigma^{2}x^{2}$ are not time-varying. I use $t$ from 0 to 1, $x$ from 0.1 to 3 and both with 100 points on grid. Matlab seems to not able to handle a matrix exceed $100000\times 100000$, so the maximum I tried so far is $15000 \times 15000$.

The test does the following things:

1. Compute $A{b}$ and $v{b}$ which is the baseline uniform grid time-varying utility case (Though the code for uniform and non-uniform are the same in my function).

2. I add 15 points both on $x$ and $t$ grid and interpolate the resulted value function on to the baseline grid. The maximum difference is $0.018$. I don't know whether this is an acceptable value, the value function ranges from 20 to 180.

Does this mean we need implicit method?

[sevhou]

Jesse: I've updated the function, with the code changed it runs with no problem, what is the problem you came across?

I'll wait till you can run the test to push the new tests to the server.

[jlperla]

@sevhou Seems to be running now. Get the latest, where I just changed one more function. Take a look at my last two commits and you will see the issue: some of the matrices were ending up dense rather than sparse.

With this, see if you can crank up the number of grid points, which should be much more doable with sparse matrices. Try to add in the grid points in the x and the t grid separately. If adding 15 points in the x grid makes 0.018 difference, that is a problem. If it does in the $t$ space, then that may make sense.

[sevhou]

@jlperla I've used sparse version of code, and now I can work on grid with size upto $1000000\times 1000000$. I've tried the new add point test, with $I=1000$ and $N=100$ (time points). Adding 15 points on x grid gives max difference on a level of $1e-4$, add 15 points on t grid gives max difference of $0.0036$, add points on both dimension gives max diff of $0.0038$.

[jlperla]

That sounds well within reasonable bounds to me (since N=100 is not all that big). But with a little more testing I think we have what we need. (I may play around with the generation of the matrix to make it a little more efficient, so hopefully the test cases will detect if I break anything).

My broader concern is convergence in general, since I suspect some sort of CFL condition might need to hold. I think it is worth doing the algebra for the implicit setup, and seeing if it makes sense. Whether we want to code it up yet is another issue. @sevhou and @stevenzhangdx : you guys can figure out who wants to play around with the tests, and who wants to try to work through the implicit algebra in the document.

[sevhou]

@jlperla We have agreed Steven will work on the implicit algebra, I'll finish the remaining test today and join Steven once I finished.

[sevhou]

I've finished up all the test I think:

1. basic test, as described before;

2. test for the case where functions not time varying. Compare result of time varying code with non-time-varying code.

3. test for $u$ being time varying and $\mu$, $\sigma^{2}$ not time varying; add point on both grid, results change within $1e-3$ (most change happens due to t-grid change). Then I checked if rais $N$ to $1000$, the max difference becomes $1e-4$.

4. test for $\mu$ being time varying only, addpoint on both grid. result change on level of $1e-4$.

5. test for both $\mu$ and $u$ being time varying, and shift most of the points on the grid, with $I=1000$, $N=1000$, the max difference between uniform grid result and shifted grid result is on level of $1e-7$.

[jlperla]

This seems to be complete. See my latest code changes for more efficient generation of the banded matrices (which still passes the test suite).