jlperla / continuous_time_methods Goto Github PK

From Ben:
Simply fix one element of the vector f, then solve the linear system (which is no longer singular) and the renormalize f afterwards. See e.g. this code snippet from http://www.princeton.edu/~moll/HACTproject/huggett_partialeq.m

AT = A';
b = zeros(2*I,1);

%need to fix one value, otherwise matrix is singular
i_fix = 1;
b(i_fix)=.1;
row = [zeros(1,i_fix-1),1,zeros(1,2*I-i_fix)];
AT(i_fix,:) = row;

%Solve linear system
gg = AT\b;
g_sum = gg'*ones(2*I,1)*da;
gg = gg./g_sum;

The question will be whether this approach can be generalized in a nonlinear setup.

See http://www.princeton.edu/~moll/HACTproject/HACT_Additional_Codes.pdf and http://www.princeton.edu/~moll/HACTproject/HACT_Numerical_Appendix.pdf

Build up to the http://www.princeton.edu/~moll/HACTproject/HJB_NGM_implicit.m solution for neoclassical growth with everything deterministic. No need to follow the explicit solution methods.

See response on https://math.stackexchange.com/questions/2466489/solving-a-large-sparse-asymmetric-linear-complementarity-problem

It seems like the discretized operator is Toeplitz as long as (1) the parameters don't change as a function of $z$; (2) the boundary conditions are `truncated' as opposed to directly putting in the reflecting barrier; and (3) grid spacing is uniform.

The truncation of the operator seems to be used in https://ecommons.cornell.edu/handle/1813/5453 and it is worth seeing if the solution is close enough. This would be unlikely to work for the stationary distribution, though.

If so, then Toeplitz solvers could be used for the sparse stationary system, and possibly even the time-varying system.

Check all of the math in the docs/operator_discretization_finite_differences.tex file. This has been separated from the older file, since I realized these notes will apply to all sorts of other examples.

A few things to check carefully:

• I added in the boundary values for both the reflection and the absorbing barrier. Go through the algebra very carefully to make sure I didn't make a mistake in them.
• You will also see that I think the formulas for (7) to (11) simplify in the boundary conditions (especially in the reflection case).
• Do a sanity check that the upwind procedure (i.e. the sign of of $\mu_1$ doesn't effect the reflection boundary values in (17) and (20). My suspicion is that it does not since this seems to make the math add up for each row adding to 1 - as it should as a proper markov chain.
• I think I modified your case of the $\mu < 0$ specialization correctly in (21) to (28), but I would check it out to be sure.

I think there are numerical instabilities being introducing for large numbers of grid points due to the division by $\Delta$ and $\Delta^2$ in the definitions.

I think I made all of the changes to the algebra document, though you should check very carefully in conjunction with #11 and then #13. After the algebra is verified, and before you do $13, you should convert the discretization code to match the algebra in the document.

Basically, the discretize_univariate_diffusion.m file just needs to match the new algebra for the A to return, and the checks in the test suite are modified to the new A.

Base this off of #8, though it could easily be done with the LCP.m based one.

Essentiallly, just do the finite differences with a fixed, but non-uniform grid. The key step will be that dx is a vector of first differences of the x vector (provided by the user).

Then when creating the X,Y,Z matrices or generating the $A$ matrix, care will need to be taken to divide by ./dx The 2nd difference operator changes as well: Care will be need to be taken if diving by dx(i) vs. dx(i+1) vs. dx(i-1) given the current setup.

See http://www.princeton.edu/~moll/HACTproject/HACT_Numerical_Appendix.pdf page 22 for most of the work done and for verification.

@FernandoKuwer Once I get this up and running, you can start adding in your tests of the DiffEqOperators.

https://github.com/chkwon/PATHSolver.jl

http://pages.cs.wisc.edu/~ferris/path/LICENSE

Following the basic pattern of #13, we want to create a set of regression tests for the simple optimal stopping problem. You will see simple_optimal_stopping_diffusion_test.m which already has a test of the code against the old HACT setup.

Some of the tests to add are:

• Some variations on the $u(x)$ function, including something convex
• A $u(x)$ function which is negative for small x
• A $S(x)$ function which is negative
• A $S(x)$ function which is so low that the agent would NEVER choose to stop.
• A $S(x)$ function which is increasing in x,
• A $S(x)$ function which is decreasing in x
• A $S(x)$ function which is non-monotonic in x, but where there is still only a single stopping point.
• A $S(x)$ function which is non-monotonic in x, but where you end up with two stopping points. You may need to mess with the $u(x)$ to get it to work.
• Solution with $\mu(x)$ always negative, always positive, and = 0
• Solution with a $\mu(x_{\min}) < 0$ and a $\mu(x_{\max}) > 0$. Then the reverse.
• Solution with $\mu(x) < 0$ and $\sigma^2 = 0$ for all $x$.
• Solution with $\mu(x) > 0$ and $\sigma^2 = 0$ for all $x$.

As always, the first part of the test is to make sure the results make sense, then the second part is to ensure it can automatically test that nothing changed.

Take a look at the stopping problem: http://www.princeton.edu/~moll/HACTproject/option_simple.pdf

Download http://www.princeton.edu/~moll/HACTproject/option_simple_LCP.m and http://www.princeton.edu/~moll/HACTproject/LCP.m

The goal is to figure out what the boundary value used for the right hand side of the finite difference scheme.. If the diffusion term was not there, this wouldn't be a problem since it uses the backward differences as upwind. But with the 2nd order term we need something else for the final point?

For the solution to the stationary_distribution_discretized_univariate.m and others, try to see if other preconditioning strategies work well for solving the LLS problem. For example, look at http://epubs.siam.org/doi/pdf/10.1137/140975358 and http://www2.cs.cas.cz/~tuma/ps/bru_mas_marin_tuma_2014.pdf for some references.

@stevenzhangdx you will see a change to the code, along with a new file run_examples.m that can be modified to drive the examples.

This will be separated from the tests subdirectory.

Want the $r$ to change smoothly.

Check for some interior set of the $x$ values (e.g i in the middle half of the statespace). FInd the
$$
\max_{1/4 I < i < 3/4 I}{v^{n+1}_i - v^n}
$$
and see how that number changes with $n$

In particular, the X matrices defined in the document are of size $(I+1)\times I$ and $(2 I + 1)\times 2I$, but the rank of them is $I$ and $2I$. This means that the overdetermined system can be solved as a linear system. The rank can be checked with an SVD or the rank function of matlab. The fact that this is really a system should be pointed out in the document prior to describing it as a LLS problem.

On the other hand, this doesn't necessarily change the solution technique for the system. For dense setups, the LLS approach can be compared to doing a pseudo-inverse or something similar in speed - but I suspect that LLS is close enough.

For sparse systems, we could investigate whether there are any sparse direct solvers that can solve rectangular systems. Iterative methods wouldn't work, but gaussian elimination should. Alternatively it is possible that sparse LLS is the best approach even if the rank is known.

When #33 is done, we can code up a basic test suite (along the lines of the current one) for the implicit system to compare it.

This is a continuation of #31

The trickiest part of this is figuring out what to do at the $t=0$ and $t=T$ steps. On the right hand side, we want this to nest solving the stationary distribution. On the left hand side, we need to remember that there is no initial condition.

After #1 is complete, you will see that the method uses an iterative step to update the value function. This works, but in principle the derivative information could provide better convergence of the model (and enable writing it in a way that can be extended).+

To do this, we first need to rewrite this as a system of equations in the vector V where a solution is given if the stacked equations are all close to 0. Unlike the existing code, which solves a linear system, this will be written as a nonlinear system of equations).

I want to update the neoclassical growth model to be as simple as possible and to test removing the dependence on the stationary solution. The task list is:

• Change all notation to match HJBE_discretization.pdf as much as possible
  • e.g. s -> gamma
  • df to Delta or whatever is the correct match)
  • min(mub,0) instead introduced mu_B_m max(muf,0)introducesmu_F_p`, etc.
• Add comment pointing to the equation number in the latex file for each line of code involving algebra. If you need to add in more equations (or want to rearrange equations) in that document, then go for it.
• We want to write the equations of the model pretending that we don't know the steady-state. To do this, we will need to play around with what is currently equation (19) in the notes.
  • First, we need to ensure that all cases of $\mu_F$ and $\mu_B$ signs are properly dealt with in the logic, as I am not sure the code/algebra is correct as stated. In particular, it should be clear what happens if both $\mu_F &gt; 0$ and $\mu_B &lt; 0$ at the same time in (19) as opposed to a footnote. Next, what happens if the drifts are exactly 0?
  • I believe that one way to write it where nothing falls through the cracks is as follows:
    $
    {v_i^n}'={v_{i,F}^n}'\bold{1}{{ \mu{i,F} \geq 0}}+{v_{i,B}^n}\bold{1}{{ \mu{i,B}\leq 0 \text{ and } \mu_{i,F} < 0 }}+{\bar{v}i^n}'\bold{1}{{ \mu_{i,F}<0<\mu_{i,B} }}\label{eq:dv-upwind}
    $
  • See if that changes the solution.
  • Next, lets see what happens if we change the behavior around the steady state to pretend we don't know it. In particular, one thing to try is whether (with a fine enough grid) we can just use $\mu_{i,F}$? In that case it would be something like
    $
    {v_i^n}'={v_{i,F}^n}'\bold{1}{{ \mu{i,F} \geq 0}}+{v_{i,B}^n}\bold{1}{{\mu{i,F} < 0 }}
    $
  • Then check if it also works for taking the average. i.e. something like
    $
    {v_i^n}'={v_{i,F}^n}'\bold{1}{{ \mu{i,F} \geq 0}}+{v_{i,B}^n}\bold{1}{{ \mu{i,B}\leq 0 \text{ and } \mu_{i,F} < 0 }}+{\frac{v_{i,F}^n + v_{i,B}^n}{2}}'\bold{1}{{ \mu{i,F}<0<\mu_{i,B} }}\label{eq:dv-upwind}
    $
  • In all cases, when you check if it "works", you should be seeing if the solution changes much, and (especially) if the correct steady state comes out naturally. Also, play around to see if it converges for different initial guesses on the value function.
  • Leave (commenting out whichever one you aren't using) code for the 3 cases in the file, so we can play with it a bit

I want to ensure that /docs/HJBE_discretization.tex has a complete set of algebra for the system of equations. I would like to write things in the notation of /docs/operator_discretization_finite_differences.tex For example, using $\mu^+$ and the $X,Y,Z$ matrices for the upwind procedure - though that will simplify due to the lack of diffusion.

After this, #2 becomes possible as the next step.

Given #16, stack with the $(\rho I - A) v = u$ equation and the system of the KFE to jointly solve with linear least squares.

The issue is that it is it is a singular matrix, and ultimately an eigen problem. However, with an overconstrained system (adding the sum to 1 constraint and the weak positivity constraints), perhaps sparse linear least squares can solve the problem? Depends on #15

With a linear operator $A$, the adjoint $A*$ can be used to find the stationary KFE. This can be an eigenvalue problem.... maybe the eigenvector associated with the unit eigenvector of $\exp(A*)$? Normalized to be positive and sum to 1?

Finish the annotation of: examples/continous_time_methods/optimal_stopping/option_simple_LCP.m

You will see a bunch of TODOs in the file. A few of them:

• Find out in both the HACTproject code and this one if it looks like it is solving geometric brownian motion, or something different. In particular, see the original code

mu_bar = -0.01;
mu = ones(I,1).*mu_bar;

sig_bar = 0.01;
sig = x.*sig_bar;
sig2 = sig.^2;

Are we right to think that this says the drift is -0.01 and the variance squared is (0.01 *x)^2 ? Why no x on the drift?

• Better explain the X, Y, and Z matrices for doing the upwind scheme in equations and math. I would like to map this into the notation of http://www.princeton.edu/~moll/HACTproject/HACT_Numerical_Appendix.pdf pages 15 and 16. In fact, those x, y, and z might be EXACTLY those with a map between s and mu If so, put in comments to explain these (and the s_F and s_B, mappings) in comments. Alternatively, perhaps the 2nd order terms are better done
• From these, explain the construction of the A matrix using spdiag better using comments, etc. Perhaps mapping to the http://www.princeton.edu/~moll/HACTproject/HACT_Numerical_Appendix.pdf setup... All of these mappings should go in the optimal stopping notes
• Better explain how the right hand boundary value is constructed with the A(N_x,N_x)= Y(N_x) + sigma_2(N_x)/(2*dx_2); A(N_x,N_x-1) = X(N_x); variation in the code.
• For the left hand boundary value: verify in your mind that it really is v(0) = 0??? and see if this is enough, or if it should be adjusted for v(x_min) = B for some $B &lt; S(x_min)$. Alternatively, the stopping may simply never bind but the dependence on the $v(0) = 0$ would be unclear.

After your modification and tests in discretize_univariate_diffusion_test.m some of the other files are broken because of the change to multiplying $A$ by $\Delta$.

Can you fix this by modifying the following functions:

• /matlab/lib/simple_HJBE_discretized_univariate.m
• /matlab/lib/simple_joint_HJBE_stationary_distribution_univariate.m
• /matlab/lib/simple_optimal_stopping_diffusion.m
• /matlab/lib/stationary_distribution_discretized_univariate.m
• /matlab/examples/optimal_stopping/LCP_example.m

The modifications will be pretty simple. In particular, mostly multiplying by $\Delta$ where currently there is an $A$ to follow the notation and equation numbers in operator_discretization_finite_differences.pdf and optimal_stopping.pdf

You will know you are successful when you can run all of the tests in the matlab folder (i.e. run_tests.m) with no errors. Also, the /matlab/examples/optimal_stopping/LCP_example.m will require running manually.

Implement the #2 function. This will involve modification of http://www.princeton.edu/~moll/HACTproject/HJB_diffusion_implicit_RBC.m and http://www.princeton.edu/~moll/HACTproject/HJB_no_uncertainty_implicit.m

A few things:

• Modifying the code, a key difference will be that do not solve the sparse system of equations for the V_stacked and then reshape to V. Instead, the "solution" of that system should be implicit in the calculation of the residual.
• Getting it to auto-differentiate may take some care, as I am not sure if reshape, spdiags, etc. directly support AD without some work. You can use loops instead as required to setup the equations.

This includes the logic on the boundary values and non-uniform grid, which you may not have run into.

Just to get started. Close once you think that you have reviewed it.

@stevenzhangdx I am starting on this, but have added you as an observer. May hand over to you for testing/etc.

While the code is written for generic upwind and downwind processes, the baseline has the same sign for all $x$.

The question is:

• what exactly are the boundary values and how do they manifest into A and how are they interpreted in the case of $\mu &gt; 0$ vs. the one we have been focusing on with $\mu &lt; 0$. What changes and why?

For example, in the $\mu &lt; 0$ the complication at the right hand side was that for a maximum of $N$ points, it was trying to reference the $N+1$ due to the 2nd derivative. The first derivative was upwind and using backward differences so didn't reverence the $N+1$ directly.

• Is it the same story with a $\mu &gt; 0$? Does it use EXACTLY the same assumptions (but now working forward? We need to write this down explicitly.

• Similar complications occur on the left hand boundary when done carefully.

• For all cases, write things up very explicitly in optimal_stopping_notes

You will see that I modified /matlab/lib/simple_HJBE_discretized_univariate.m (in 76e9373 ) 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).

Once #5, #6, and #7 are complete, we can swap out the LCP function used in the baseline example for a high-performance commerical version. This is unlikely to help with this exact scenario since it is too simple, but may be necessary for larger problems with many more dimensions.

With tomab, see https://tomopt.com/docs/quickguide/quickguide027.php as a starting point. Knitro is likely to be the best solution.
It is essential that you verify it is setup (and solving) as a sparse problem. May require flipping through documents to look for possible settings that ensure sparseness, and setting sparseness patterns directly (or indirectly)?

Worst case, any quadratic solver with linear constraints can do LCP, since there is a mapping of the LCP to the KKT conditions of a quadratic problem. See https://en.wikipedia.org/wiki/Linear_complementarity_problem#Convex_quadratic-minimization:_Minimum_conditions for an example of the mapping.

Relies on #40 for the basic package.

You will need to get your basic julia workflow working with Juno/etc. Take a look at https://github.com/econtoolkit/julia

From the univariate diffusion, get the discretized operator with the reflecting boundaries, and compare to the matlab versions. We will add automatic tests later.

You will see that this is a split file which takes the discretization of the operator as given. Check the algebra against http://www.princeton.edu/~moll/HACTproject/option_simple.pdf and make sure that I added the slack variable correctly in my addition, as we may need to use that specification for some approaches.

Talk to Chris R about this as required to get the full DiffEqOperators setup. In the background, worth watching his tutorial to get a sense of how this stuff works: https://www.youtube.com/watch?v=tx8DRc7_c9I&amp=&t=2s

For our own simple testing setup, I will leave things much more open.

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.

Hi guys: I have set you both as assignees, and you can discuss splitting up any tasks.

I think I have a starting point for the algebra for discretizing the PDE. Get the latest version of the operator_discretization_finite_differences.tex and compile it. You will see section 2 now has code on discretizing an operator with time dependence. Hopefully it is clear.

Some of the tasks here are:

• Verify that the algebra on (69) to (77) is correct. This is basically just doing the same thing we did for the stationary setup, but with the potentially time varying drift, etc. I think everything should be at $n$.
• Fill I suspect that (80) and (81) are relatively close to correct, but there is nothing done for the corners (which now include both $\underline{z}$ and $\bar{z}$ along with $t=0$ and $t=T$. Hopefully we can use the BC in (66). For the $t=0$, if we use an explicit time step, I am not sure we need anything special? For implicit time steps things get more confusing.
• See if (82) is correct, and fix it otherwise. This was purely a guess, but I suspect the correct answer is something close to this.
• If it is helpful, you can do a uniform-grid version of all of these, but my guess is that it won't make life sufficiently easier to be worth it.

Can we get analytical gradients with: http://tomsym.com/index.html

You will see a starting point in continuous_time_methods\matlab\tests\discretize_univariate_diffusion_test.m
The goal is to fully test the continuous_time_methods\matlab\lib\discretize_univariate_diffusion.m script.

To come up with the comprehensive test suite, you will want to go through all of the key permutations of parameters. The idea is that everytime we make changes to the discretize_univariate_diffusion.m file, we run the regression test to make sure nothing changed.

Beyond just ensuring that nothing changed, you should also make sure that the results match the theory in the continuous_time_methods\docs\optimal_stopping.pdf notes. This is sometimes best done with simple versions (e.g. 5 nodes) where you can verify by hand that things make sense.

Some of the parameter variations to check:

• x_min < 0
• Changes in the scale of x_min and x_max give the correct answer for a given $I$
• $\mu(x) = 0$ for all $\mu$
• $\mu(x) = \mu &gt; 0$ to make sure that the upwind approach makes sense. Special care to think through the corners!
• $\mu(x) = \mu &lt; 0$ to make sure that the upwind approach makes sense. Special care at the corners.
• $\mu(x_min) &lt; 0$ and $\mu(x_max) &gt; 0$ monotonically increasing.
• $\mu(x_min) &gt; 0$ and $\mu(x_max) &lt; 0$ monotonically decreasing.
• $\mu(x_min) &lt; 0$ and $\mu(x_max) &lt; 0$ concave so that $\mu(x) &gt; 0$ at some point. Make sure that the finite-differences locally to that sign change makes sense and matches what you would calculate by hand!
• $\sigma^2 = 0$ everywhere.
• $\sigma^2 = 0$ and varying in a few places.

Given #37, you should fill in the continuous_time_methods\matlab\tests\time_varying_optimal_stopping_diffusion_test.m tests with important tests to make sure you think the method is correct, and delivers sensible results. I threw in a few starting points, but they are not really verified.

• Check that if we have constant $\mu,\sigma^2,S,u$ but multiple time periods, that the value functions and stopping point are the same.
• Have the S function increasing over time, and make sure that both (a) the $x$ at which stopping occurs is increasing, but also that the point at which stopping occurs is strictly greater than it is for a constant S at that point (i.e., agents would wait around for the S to grow and stop less)
• Do the same with a decreasing S
• Do the same with a $u(t,x)$ that is time varying. Ensure it moves large enough that we get more than trivial changes in the soluion.
• Try to have a time-varying $\mu(t,x)$ and make sure that the solution is sensible from an agents perspective.

Come up with a test case to ensure no bugs were introduced into the new code.

Save a vector of out put from original HACT code, put it in folder, and then automatically test at the bottom of the new file if it is the same.

@stevenzhangdx adding as observer. This depends on #28

First, go through the optimal_stopping.pdf section 1 to make sure that you understand the stopping problem and how it maps to a LCP.

Then, go through the new section 2 to make sure I didn't introduce a mistake. It is very simple, as it relies on our discretization of the time-varying diffusion process (which we have hopefully fully tested by this point).

Finally, you will see in continuous_time_methods\matlab\lib\optimal_stopping_diffusion.m an implementation of this method, which relies on the discretization of the linear operators. I wrote up a starting point for the tests in time_varying_optimal_stopping_diffusion_test.m. Make sure that you think this implements what it says it does, and that the code matches the description. At this point, no need to test anything besides the yuval algorithm