Improve speed of ConsGenIncProcessModel solvers

[sbenthall]

The ConsGenIncProcessModel solvers, with default values, are currently the longest running automated tests.

Maybe there is a way to improve the performance of these solvers.

[mnwhite]

Why not just have them only solve for a few periods? Just set cycles=1 or 5 or something.

On Thu, Mar 12, 2020 at 5:46 PM Sebastian Benthall notifications@github.com wrote:

The ConsGenIncProcessModel solvers, with default values, are currently the longest running automated tests.

Maybe there is a way to improve the performance of these solvers.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/econ-ark/HARK/issues/567, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADKRAFMEFDIRJ5XTDLIGWE3RHFJ2RANCNFSM4LGV5CGQ .

[sbenthall]

Hmmm. Ok, I'm trying that,

So, with cycles=1, it looks like the solver does not make much progress on converging on a solution.

With cycles=5, I get errors in the updatepLvlGrid method, because of an assertion the cycles equal either 0 or 1.

https://github.com/econ-ark/HARK/blob/d1a8eb90ebe34a93801ead8655b5e29ef9d1e3eb/HARK/ConsumptionSaving/ConsGenIncProcessModel.py#L1142

[sbenthall]

I guess I'll go with cycles=1.

I wonder if you could help me understand what's going on here though. Why is does cycles=1 mean that the model is a lifecycle model?

[mnwhite]

Sure. The cycles attribute indicates how many times the sequence of periods described by the time-varying parameters will "happen" to the agent, with 0 standing in for infinity. For any value of cycles other than 0, it's a finite horizon or lifecycle model-- there is some definite end point that the agent cannot exist past. Basically every model that everyone ever works with would be specified as either cycles=0 or cycles=1; having cycles >= 2 isn't completely useless, but it's best used for teaching (or sometimes debugging). For example, you can quickly plot how an infinite horizon solution converges by setting cycles=50 and solving a problem with a 1 period cycle (the policy function visually converges in about 50 periods).

The default dictionary for ConsGenIncProcess model specifies a cycle that is T_cycles=1 periods long. If you set cycles=0, the model solves the problem of "what should I do if I will live out this kind of period over and over again?". If you set cycles to 1 with a one period cycle, you are asking it to solve the problem of "what should I do if I live this period, and then it's the terminal period"? The solution isn't trying to converge, you're just asking it to solve a problem with one non-terminal period and one (trivial) terminal period. It's a lifecycle model for a mayfly.

Setting cycles=5 works fine for most models in HARK, but ConsGenIncProcess is a special case because of the pLvlGrid. This attribute is constructed in one of the methods called by update; it represents the grid of permanent income levels at which to solve for optimal consumption, with the consumption function interpolated in between. It's constructed by simulation, literally drawing out the population of pLvl at each time index of the cycle and then grabbing percentiles as given in the primitive attribute pLvlPctiles. In the default dictionary, it captures the middle 99.8% of the distribution. If cycles=1, this just means looping over each period of the lifecycle and finding the distribution of permanent income at each age. If cycles=0, the code assumes that 1000 periods is a "long time" and simulates the population for that long, then uses the distribution of pLvl in that fake-ergodic distribution at each t of the cycle.

But if cycles >= 2, this method won't work. The distribution of pLvl is different on successive cycles; there is not a one-to-one mapping between t_age and t_cycle, nor is there a meaningful "long run" distribution. There are ways around this, but I didn't take a stand on how we should handle this non-standard case when I wrote that method. As is, if you really wanted to have a sequence of periods (or just one period) repeat a finite number of times, you would have to write the parameter dictionary to have that loop in every time-varying primitive parameter, and then set cycles=1. For example, PermGroFac would change from being [1.02] to 5*[1.02].

On Fri, Mar 13, 2020 at 12:35 PM Sebastian Benthall < notifications@github.com> wrote:

I guess I'll go with cycles=1.

I wonder if you could help me understand what's going on here though. Why is does cycles=1 mean that the model is a lifecycle model?

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/econ-ark/HARK/issues/567#issuecomment-598810289, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADKRAFIIXBZB2H2UXZK7IE3RHJOFVANCNFSM4LGV5CGQ .

[sbenthall]

Thanks for explaining all that.

So if I understand correctly...

• T_cycles is the number of periods in a cycle, a property of the model
• cycles is the number (probably infinite or 1) of times the agent lives through the cycle, for the purpose of the solver
• T_sim is the number of periods to be simulated when simulate() is called

And ConsGenIncProcess doesn't use dynamic programming in the solver in the same way as the other models, but uses simulation internally, with a 1000 period limit.

I've made some changes to #527 that set cycles=1 for the tests, but these tests at this point just make sure the code doesn't break; it doesn't get at the internals at all.

But maybe there is a way to improve the speed of the ConsGenIncProcessModel solver by introducing a threshold that would interrupt the solver's internal simulation if the model converges before 1000 periods. It sounds like 50 is sufficient to prove the concept?

[mnwhite]

Not quite. Both T_cycle and cycles are instance-level attributes. The model that an AgentType subclass solves is agnostic to T_cycle and cycles; the specific instance of that model (a "problem") fills in those time-defining attributes.

T_sim as an attribute represents the maximum number of periods you would ever be interested in simulating for this instance. It's used to initialize the arrays that will track idiosyncratic variable histories of the attributes named in track_vars. When simulate() is called, it will simulate for T_sim period by default, but you can also specify an argument T <= T_sim for it to simulate for. This is useful for doing things like "simulate 100 periods, then screw with the state of the world, then simulate 300 more periods".

ConsGenIncProcess uses backward iteration methods to solve the model, just like all our other solvers. It simply uses simulation to construct the grid the problem will be solved on, in the pre-solution step. The 1000 periods thing is just an approximation to the "long run" distribution of permanent income. I.e. if a model period is a year, it assumes that after 1000 years, the population distribution of pLvl is roughly distributed the same as it would after 10,000 years or 10B years.

"Making sure the code doesn't break" is a really, really important test. Having the solver try a 1 or 2 period problem and then check to see if a few points of the policy function match the solution from the prior version of the code is pretty much the best we can do with our model tests.

Your last paragraph is incorrectly premised for the reasons above.

On Fri, Mar 13, 2020 at 2:33 PM Sebastian Benthall notifications@github.com wrote:

Thanks for explaining all that.

So if I understand correctly...

• T_cycles is the number of periods in a cycle, a property of the model
• cycles is the number (probably infinite or 1) of times the agent lives through the cycle, for the purpose of the solver
• T_sim is the number of periods to be simulated when simulate() is called

And ConsGenIncProcess doesn't use dynamic programming in the solver in the same way as the other models, but uses simulation internally, with a 1000 period limit.

I've made some changes to #527 https://github.com/econ-ark/HARK/pull/527 that set cycles=1 for the tests, but these tests at this point just make sure the code doesn't break; it doesn't get at the internals at all.

But maybe there is a way to improve the speed of the ConsGenIncProcessModel solver by introducing a threshold that would interrupt the solver's internal simulation if the model converges before 1000 periods. It sounds like 50 is sufficient to prove the concept?

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/econ-ark/HARK/issues/567#issuecomment-598858948, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADKRAFJDEP5TJGGWNCTNBLTRHJ36XANCNFSM4LGV5CGQ .

[mnwhite]

Also, I just noticed that our documentation notebooks are in the /examples folder (at least on the branch I'm working on), and the actual examples have been deleted and replaced with bizarro notebook output. Is this intended? I thought the whole point of the /examples directory is to have small, minimal .py files that implement our models?

On Fri, Mar 13, 2020 at 2:48 PM Matthew White mnwhite@gmail.com wrote:

Not quite. Both T_cycle and cycles are instance-level attributes. The model that an AgentType subclass solves is agnostic to T_cycle and cycles; the specific instance of that model (a "problem") fills in those time-defining attributes.

T_sim as an attribute represents the maximum number of periods you would ever be interested in simulating for this instance. It's used to initialize the arrays that will track idiosyncratic variable histories of the attributes named in track_vars. When simulate() is called, it will simulate for T_sim period by default, but you can also specify an argument T <= T_sim for it to simulate for. This is useful for doing things like "simulate 100 periods, then screw with the state of the world, then simulate 300 more periods".

ConsGenIncProcess uses backward iteration methods to solve the model, just like all our other solvers. It simply uses simulation to construct the grid the problem will be solved on, in the pre-solution step. The 1000 periods thing is just an approximation to the "long run" distribution of permanent income. I.e. if a model period is a year, it assumes that after 1000 years, the population distribution of pLvl is roughly distributed the same as it would after 10,000 years or 10B years.

"Making sure the code doesn't break" is a really, really important test. Having the solver try a 1 or 2 period problem and then check to see if a few points of the policy function match the solution from the prior version of the code is pretty much the best we can do with our model tests.

Your last paragraph is incorrectly premised for the reasons above.

On Fri, Mar 13, 2020 at 2:33 PM Sebastian Benthall < notifications@github.com> wrote:

Thanks for explaining all that.

So if I understand correctly...

• T_cycles is the number of periods in a cycle, a property of the model
• cycles is the number (probably infinite or 1) of times the agent lives through the cycle, for the purpose of the solver
• T_sim is the number of periods to be simulated when simulate() is called

And ConsGenIncProcess doesn't use dynamic programming in the solver in the same way as the other models, but uses simulation internally, with a 1000 period limit.

I've made some changes to #527 https://github.com/econ-ark/HARK/pull/527 that set cycles=1 for the tests, but these tests at this point just make sure the code doesn't break; it doesn't get at the internals at all.

But maybe there is a way to improve the speed of the ConsGenIncProcessModel solver by introducing a threshold that would interrupt the solver's internal simulation if the model converges before 1000 periods. It sounds like 50 is sufficient to prove the concept?

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/econ-ark/HARK/issues/567#issuecomment-598858948, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADKRAFJDEP5TJGGWNCTNBLTRHJ36XANCNFSM4LGV5CGQ .

[mnwhite]

Nevermind, I found it one folder down.

Open up examples/ConsumptionSaving/example_ConsGenIncProcessModel.py ... At lines 38 and 130 you can add the argument True to the call to solve() if you want to see it run cycle by cycle. These examples only take about 6-8 seconds to solve; this doesn't seem like a long time to me.

Also, it looks like process_time is giving weird numbers. Can we use time.time instead?

On Fri, Mar 13, 2020 at 2:56 PM Matthew White mnwhite@gmail.com wrote:

Also, I just noticed that our documentation notebooks are in the /examples folder (at least on the branch I'm working on), and the actual examples have been deleted and replaced with bizarro notebook output. Is this intended? I thought the whole point of the /examples directory is to have small, minimal .py files that implement our models?

On Fri, Mar 13, 2020 at 2:48 PM Matthew White mnwhite@gmail.com wrote:

Not quite. Both T_cycle and cycles are instance-level attributes. The model that an AgentType subclass solves is agnostic to T_cycle and cycles; the specific instance of that model (a "problem") fills in those time-defining attributes.

T_sim as an attribute represents the maximum number of periods you would ever be interested in simulating for this instance. It's used to initialize the arrays that will track idiosyncratic variable histories of the attributes named in track_vars. When simulate() is called, it will simulate for T_sim period by default, but you can also specify an argument T <= T_sim for it to simulate for. This is useful for doing things like "simulate 100 periods, then screw with the state of the world, then simulate 300 more periods".

ConsGenIncProcess uses backward iteration methods to solve the model, just like all our other solvers. It simply uses simulation to construct the grid the problem will be solved on, in the pre-solution step. The 1000 periods thing is just an approximation to the "long run" distribution of permanent income. I.e. if a model period is a year, it assumes that after 1000 years, the population distribution of pLvl is roughly distributed the same as it would after 10,000 years or 10B years.

"Making sure the code doesn't break" is a really, really important test. Having the solver try a 1 or 2 period problem and then check to see if a few points of the policy function match the solution from the prior version of the code is pretty much the best we can do with our model tests.

Your last paragraph is incorrectly premised for the reasons above.

On Fri, Mar 13, 2020 at 2:33 PM Sebastian Benthall < notifications@github.com> wrote:

Thanks for explaining all that.

So if I understand correctly...

• T_cycles is the number of periods in a cycle, a property of the model
• cycles is the number (probably infinite or 1) of times the agent lives through the cycle, for the purpose of the solver
• T_sim is the number of periods to be simulated when simulate() is called

And ConsGenIncProcess doesn't use dynamic programming in the solver in the same way as the other models, but uses simulation internally, with a 1000 period limit.

I've made some changes to #527 https://github.com/econ-ark/HARK/pull/527 that set cycles=1 for the tests, but these tests at this point just make sure the code doesn't break; it doesn't get at the internals at all.

But maybe there is a way to improve the speed of the ConsGenIncProcessModel solver by introducing a threshold that would interrupt the solver's internal simulation if the model converges before 1000 periods. It sounds like 50 is sufficient to prove the concept?

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/econ-ark/HARK/issues/567#issuecomment-598858948, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADKRAFJDEP5TJGGWNCTNBLTRHJ36XANCNFSM4LGV5CGQ .

[llorracc]

@sbenthall please confer with @MridulS who wrote the excellent configurator.py tool for trying out various configurations of parameter values for the CGMPortfolio REMARK. This was intended as a prototype for a more general purpose tool that I hope we can construct that will define a standard procedure for configuring and running our models. That way, there could be a special configuration for DemARK testing that, for example, had just a 5 or 10 period configuration for ConsGenIncProcessModel instead of its default. We could have special "error-testing" configurations for all DemARKs (and for that matter REMARKs) that would be constructed to run fast while still doing at least a bit of a workout of the code.

[sbenthall]

Keeping to the original question:

ConsGenIncProcess uses backward iteration methods to solve the model, just like all our other solvers. It simply uses simulation to construct the grid the problem will be solved on, in the pre-solution step. The 1000 periods thing is just an approximation to the "long run" distribution of permanent income. I.e. if a model period is a year, it assumes that after 1000 years, the population distribution of pLvl is roughly distributed the same as it would after 10,000 years or 10B years.

So, suppose one wanted to write a test for this functionality, but did not want to wait for it to go through 1000 steps.

Maybe there is or could be a way to limit that number of steps to 50, for the purpose of testing.

[mnwhite]

As written, no. The idea that 1000 periods is a "long time" in model terms is hard coded. It's not pulled from an attribute.

That could be changed with 2 lines, but there's not really a point. That construction step takes 1 second or so, probably much less. So there isn't really a point in speeding it up for testing purposes. It's literally just drawing 10,000,000 or so lognormal shocks, executing 10,000,000 multiplications, sorting 10,000 numbers, and pulling a few percentiles. This is not a burdensome step.

On Fri, Mar 13, 2020, 7:19 PM Sebastian Benthall notifications@github.com wrote:

Keeping to the original question:

ConsGenIncProcess uses backward iteration methods to solve the model, just like all our other solvers. It simply uses simulation to construct the grid the problem will be solved on, in the pre-solution step. The 1000 periods thing is just an approximation to the "long run" distribution of permanent income. I.e. if a model period is a year, it assumes that after 1000 years, the population distribution of pLvl is roughly distributed the same as it would after 10,000 years or 10B years.

So, suppose one wanted to write a test for this functionality, but did not want to wait for it to go through 1000 steps.

Maybe there is or could be a way to limit that number of steps to 50, for the purpose of testing.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/econ-ark/HARK/issues/567#issuecomment-598969671, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADKRAFK4DM7LENYREOCWILTRHK5PHANCNFSM4LGV5CGQ .

[sbenthall]

Ah, I was under the impression that this simulation step was the performance bottleneck. I may be mistaken.

The next thing I'll try is performance profiling tools to see what's slowing the solver down.

[mnwhite]

How long is the solver taking for you? It should only take a few seconds if it's running the default dictionary.

Any profiler will tell you that the most expensive line is whichever one evaluates vPfuncNext. That's where literally 95% of any computation will be, as it's evaluating a linear interpolation over 1D interpolators for about 485624 realizations of mNrmNext (if I remember gridsizes in the default dictionary correctly).

On Fri, Mar 13, 2020, 9:09 PM Sebastian Benthall notifications@github.com wrote:

Ah, I was under the impression that this simulation step was the performance bottleneck. I may be mistaken.

The next thing I'll try is performance profiling tools to see what's slowing the solver down.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/econ-ark/HARK/issues/567#issuecomment-598990572, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADKRAFPPCKODVDTQDHV4ULDRHLKLTANCNFSM4LGV5CGQ .

[llorracc]

Maybe there is or could be a way to limit that number of steps to 50, for the purpose of testing.

Another use case for the comparator-type tool that @MridulS is working on.

So, suppose one wanted to write a test for this functionality, but did not want to wait for it to go through 1000 steps.

Actually, there are a number of potential tests. These notes derive a number of propositions about the long run distribution of a process like

$p{t+1} = \gamma + \phi p{t} + \epsilon_{t+1}$

the most important of which is probably:

$\sigma^{2}{p} = \left(\frac{\phi}{1-\phi^{2}}\right) \sigma^{2}{\epsilon}$

and another of which is that $p$ itself should be distributed according to

$p \sim N\left(\left(\frac{\gamma}{1-\phi}\right),\sigma^{2}_{p}\right)$

Actually, it might be better simply to construct this distribution at the outset and then see whether after some number of periods of simulation (50?) the distribution's mean and variance are "close enough" to the analytical values (where "close enough" would be, say, only one chance in 10000 that the mean of the resulting distribution would be farther away than $x$ from its analytical value. If the number of periods is $n$ then the process of simulating $n$ periods is to add $n$ normally distributed shocks to the process, and then to do a test of whether after those shocks the data look like they are still normally distributed according to the original, supposedly stationary, distribution. (There are standard tests for whether a distribution is normal with a given variance -- I think a $\chi^{2}$ test.)

[mnwhite]

We can put in tests like that, but they shouldn't go in ConsGenIncProcessModel. The point of this model is that (log) permanent income can follow any process, not just p_{t+1} = \gamma + \phit + \epsilon{t+1} (random walk with drift; replacement/death events also occur). That special process is covered by ConsIndShockModel, and it allows us to normalize out permanent income when solving the model. And as Chris notes, he has derived closed form solutions for the long run distribution of p_t for that special process. If anywhere, these kind of tests should be built from that model.

ConsGenIncProcess exists to handle the generic case: P_{t+1} = SomeFunction(Pt) * \psi{t+1}. Note the change in notation from lower (logs) to capital (levels) for permanent income, and the switch from \epsilon to \psi. If we are completely ignorant about SomeFunction, then there is no way to know what the long run distribution of P_t or p_t should be, and thus no tests we should run. There are special cases where this can be derived, like an AR1 in logs, but there's no generic result.

On Sat, Mar 14, 2020 at 7:59 AM Christopher Llorracc Carroll < notifications@github.com> wrote:

Maybe there is or could be a way to limit that number of steps to 50, for the purpose of testing.

Another use case for the comparator-type tool that @MridulS https://github.com/MridulS is working on.

So, suppose one wanted to write a test for this functionality, but did not want to wait for it to go through 1000 steps.

Actually, there are a number of potential tests. These notes http://halweb.uc3m.es/esp/Personal/personas/amalonso/esp/TSAtema4.pdf derive a number of propositions about the long run distribution of a process like

$p{t+1} = \gamma + \phi p{t} + \epsilon_{t+1}$

the most important of which is probably:

$\sigma^{2}{p} = \left(\frac{\phi}{1-\phi^{2}}\right) \sigma^{2} {\epsilon}$

and another of which is that $p$ itself should be distributed according to

$p \sim N\left(\left(\frac{\gamma}{1-\phi}\right),\sigma^{2}_{p}\right)$

Actually, it might be better simply to construct this distribution at the outset and then see whether after some number of periods of simulation (50?) the distribution's mean and variance are "close enough" to the analytical values (where "close enough" would be, say, only one chance in 10000 that the mean of the resulting distribution would be farther away than $x$ from its analytical value. If the number of periods is $n$ then the process of simulating $n$ periods is to add $n$ normally distributed shocks to the process, and then to do a test of whether after those shocks the data look like they are still normally distributed according to the original, supposedly stationary, distribution. (There are standard tests for whether a distribution is normal with a given variance -- I think a $\chi^{2}$ test.)

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/econ-ark/HARK/issues/567#issuecomment-599047808, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADKRAFKVX2HGIFYTKQUZPR3RHNWRVANCNFSM4LGV5CGQ .

[sbenthall]

How long is the solver taking for you? It should only take a few seconds if it's running the default dictionary.

Running this test file: https://github.com/econ-ark/HARK/blob/master/HARK/ConsumptionSaving/tests/test_ConsGenIncProcessModel.py

is taking 62.7 seconds on my machine.

[mnwhite]

It looks to me like those tests are equivalent to what's being run in example_ConsGenIncProcess, but with a shorter simulation period. Those two examples take 9.2 seconds and 9.5 seconds to solve on my computer, and the simulation should take about half a second (it takes 10 seconds to simulate 500 periods).

For all of our solvers, I recommend against testing them by solving an infinite horizon problem. We should set cycles to something like 2 or 5, and target our tests off of that small-ish model. These are tests of whether the solver produces the same result after a code change.

On Sat, Mar 14, 2020 at 9:41 AM Sebastian Benthall notifications@github.com wrote:

How long is the solver taking for you? It should only take a few seconds if it's running the default dictionary.

Running this test file:

https://github.com/econ-ark/HARK/blob/master/HARK/ConsumptionSaving/tests/test_ConsGenIncProcessModel.py

is taking 62.7 seconds on my machine.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/econ-ark/HARK/issues/567#issuecomment-599062534, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADKRAFJXURZ67H6LBT3GRLDRHOCQJANCNFSM4LGV5CGQ .

[llorracc]

Because GitHub markdown does not display math, and raw LaTeX code looks terrible, @mnwhite made a guess about what I was probably saying, but it's not what I was actually saying. My math was for the specialization of the general income process that Matt mentions at the end of his comment: log p{t+1} = \gamma + \phi \log p{t} + \epsilon_{t+1}. This would be a good default for simulating income dynamics because there is a closed form solution for the steady state distribution (per my math -- copy and paste into a markdown editor that can handle LaTeX and then it becomes readable!).

And as Chris notes, he has derived closed form solutions for the long run distribution of p_t for that special process.

Actually, for the special process with purely permanent shocks and a positive probability of death, the formula in cstwMPC is for the variance of the square of P; I don't have a formula for the steady state distribution (which is unbounded above but not lognormal). But it would be easy enough to test whether the simulation produces something close to the derived variance. (I did that when I originally derived the formula to make sure my math was right, and confirmed the result).

[llorracc]

Oops, finger slipped and hit "close and comment" rather than just "comment". Reopening.