Breaking-up dr1.m

Currently the function dr1.m contains most of the code for computing first, second and third order approximation of stochastic model. Current shortcomings:

• This code is messy
• It contains parts that could be advantageously factorized so as to be computed only once per model and not repeatedly every time that one calls the solution.
• Since we limited the model to have leads and lags on one period only and added necessary auxilary variables, some part of the code can be greatly simplified.

Branch dr1break

• dr1.m is replaced by stochastic_solvers.m, but the parent function resol.m keeps calling dr_block.m separately.

• stochastic_solvers.m calls in turn

  • dyn_first_order_solver.m

  • AIM_first_order_solver.m

  • dyn_second_order_solver.m

  • k_order_pert.m

  • dyn_risky_steadystate_solver.m

• stochastic_solvers.m performs the following inside the function itself:

  • obtains first and second derivatives from <modfile>_dynamic() or bytecode()

  • solves for first order approximation of purely backward for purely forward directly
  • transforms the solution for log-linear option
  • computes the first-order coefficients of deterministic shocks

Still to be done to meet goals described below

• 2nd and 3rd order approximation of purely backward and purely forward models (making separate functions for purely backward and purely forward models)
• factorizing code for description of state space

• integrating dr_block.m in the new organization of code (separating evaluation of the Jacobian blocks and computation of the solution)

• implementing better storage of results

Proposed re-organization of the code

Factorized code

• computation of the elements of the state space and the number of different types of variables (static, forward, predetermined, both). This code should return list of variable indices.

Solution code

• replace dr1.m by stochastic_solvers.

• if k_order algorithm || 3rd order

  • k_order_perturbation

• else if 1st order

  • first_order_solvers

    • if purely backward
      • inversion of Jacobian matrix
    • if purely forward
      • algorithm to be developed (not urgent)
    • general case
      • Generalized Schur decomposition (current algorithm by default)
        • will use Sims qz routines
      • Generalized Schur decomposition for several model blocks
      • Generalized Schur decomposition (DLL version)
        • will be new default
      • Generalized Schur decomposition (DLL version) for several blocks
      • AIM algorithm
      • Partial Information
    • solution for deterministic exogenous variables
• else if 2nd order

  • second order solvers

    • current algorithm
    • approximation around the stochastic steady state (to be developed by Michel)
• Create a separate function for computing the eigenvalues of a model using only the Generalized Schur decomposition algorithm (not necessarily a separate function after all)

Storing results

• The code must keep storing results in the current format to keep compatibility with existing Dynare code and user code
• However, the current ordering of variables in stored matrices was a very bad idea, inefficient and source of many confusion. We should also save

  • matrices of derivatives of the solution function (declared variables x state variables). The variables should be ordered in declaration order (i.e. as in the VAR statement). I don't think that we need the solution function for the auxiliary variables for forward expression, but I may be wrong.

  • matrices of derivatives of the solution function for state variables (state variables x state variables). Also in variable declaration order. This is duplicate saving, but I think it is more efficient than having to reform these matrices from the previous ones over and over again.
• Systematize order_var in internal storage of matrix (let think more)

Open questions

• Shall we reorder the Jacobian to permit block copy in forming D and E