+
+| Revision |
+ Version |
+ Date |
+ Description of changes |
+
+
+
+|
+ | 1.3.7
+ | 2008/01/15
+ |
+
+ |
| | | Corrected a serious bug in centralizing a
+decision rule. This bug implies that all results based on simulations
+of the decision rule were wrong. However results based on stochastic
+fix points were correct. Thanks to Wouter J. den Haan and Joris de Wind!
+
+ |
| | | Added options --centralize and --no-centralize.
+
+ |
| | | Corrected an error of a wrong
+variance-covariance matrix in real-time simulations (thanks to Pawel
+Zabzcyk).
+
+ |
| | | Corrected a bug of integer overflow in refined
+faa Di Bruno formula if one of refinements is empty. This bug appeared
+when solving models without forward looking variables.
+
+ |
| | | Corrected a bug in the Sylvester equation
+formerly working only for models with forward looking variables.
+
+ |
| | | Corrected a bug in global check printout.
+
+ |
| | | Added generating a dump file.
+
+ |
| | | Fixed a bug of forgetting repeated assignments
+(for example in parameter settings and initval).
+
+ |
| | | Added a diff operator to the parser.
+
+ |
+| 1539
+ | 1.3.6
+ | 2008/01/03
+ |
+
+ |
| | | Corrected a bug of segmentation faults for long
+names and path names.
+
+ |
| | | Changed a way how random numbers are
+generated. Dynare++ uses a separate instance of Mersenne twister for
+each simulation, this corrects a flaw of additional randomness caused
+by operating system scheduler. This also corrects a strange behaviour
+of random generator on Windows, where each simulation was getting the
+same sequence of random numbers.
+
+ |
| | | Added calculation of conditional distributions
+controlled by --condper and --condsim.
+
+ |
| | | Dropped creating unfoled version of decision
+rule at the end. This might consume a lot of memory. However,
+simulations might be slower for some models.
+
+ |
+| 1368
+ | 1.3.5
+ | 2007/07/11
+ |
+
+ |
| | | Corrected a bug of useless storing all derivative
+indices in a parser. This consumed a lot of memory for large models.
+
+ |
| | | Added an option --ss-tol controlling a
+tolerance used for convergence of a non-linear solver.
+
+ |
| | | Corrected buggy interaction of optimal policy
+and forward looking variables with more than one period.
+
+ |
| | | Variance matrices can be positive
+semidefinite. This corrects a bug of throwing an error if estimating
+approximation errors on ellipse of the state space with a
+deterministic variable.
+
+ |
| | | Implemented simulations with statistics
+calculated in real-time. Options --rtsim and --rtper.
+
+ |
+| 1282
+ | 1.3.4
+ | 2007/05/15
+ |
+
+ |
| | | Corrected a bug of wrong representation of NaN in generated M-files.
+
+ |
| | | Corrected a bug of occassionaly wrong evaluation of higher order derivatives of integer powers.
+
+ |
| | | Implemented automatic handling of terms involving multiple leads.
+
+ |
| | | Corrected a bug in the numerical integration, i.e. checking of the precision of the solution.
+
+ |
+| 1090
+ | 1.3.3
+ | 2006/11/20
+ |
+
+ |
| | | Corrected a bug of non-registering an auxiliary variable in initval assignments.
+
+ |
+| 988
+ | 1.3.2
+ | 2006/10/11
+ |
+
+ |
| | | Corrected a few not-serious bugs: segfault on
+some exception, error in parsing large files, error in parsing
+matrices with comments, a bug in dynare_simul.m
+
+ |
| | | Added posibility to specify a list of shocks for
+which IRFs are calculated
+
+ |
| | | Added --order command line switch
+
+ |
| | | Added writing two Matlab files for steady state
+calcs
+
+ |
| | | Implemented optimal policy using keyword
+planner_objective and planner_discount
+
+ |
| | | Implemented an R interface to Dynare++ algorithms
+(Tamas Papp)
+
+ |
| | | Highlevel code reengineered to allow for
+different model inputs
+
+ |
+| 799
+ | 1.3.1
+ | 2006/06/13
+ |
+
+ |
| | | Corrected few bugs: in error functions, in linear algebra module.
+
+ |
| | | Updated dynare_simul.
+
+ |
| | | Updated the tutorial.
+
+ |
| | | Corrected an error in summing up tensors where
+setting up the decision rule derivatives. Thanks to Michel
+Juillard. The previous version was making deterministic effects of
+future volatility smaller than they should be.
+
+ |
+| 766
+ | 1.3.0
+ | 2006/05/22
+ |
+
+ |
| | | The non-linear solver replaced with a new one.
+
+ |
| | | The parser and derivator replaced with a new
+code. Now it is possible to put expressions in parameters and initval
+sections.
+
+ |
+| 752
+ | 1.2.2
+ | 2006/05/22
+ |
+
+ |
| | | Added an option triggering/suppressing IRF calcs..
+
+ |
| | | Newton algortihm is now used for fix-point calculations.
+
+ |
| | | Vertical narrowing of tensors in Faa Di Bruno
+formula to avoid multiplication with zeros..
+
+ |
+| 436
+ | 1.2.1
+ | 2005/08/17
+ |
+
+ |
| | | Faa Di Bruno for sparse matrices optimized. The
+implementation now accommodates vertical refinement of function stack
+in order to fit a corresponding slice to available memory. In
+addition, zero slices are identified. For some problems, this implies
+significant speedup.
+
+ |
| | | Analytic derivator speedup.
+
+ |
| | | Corrected a bug in the threading code. The bug
+stayed concealed in Linux 2.4.* kernels, and exhibited in Linux 2.6.*,
+which has a different scheduling. This correction also allows using
+detached threads on Windows.
+
+ |
+| 410
+ | 1.2
+ | 2005/07/29
+ |
+
+ |
| | | Added Dynare++ tutorial.
+
+ |
| | | Changed and enriched contents of MAT-4 output
+file.
+
+ |
| | | Corrected a bug of wrong variable indexation
+resulting in an exception. The error occurred if a variable appeared
+at time t-1 or t+1 and not at t.
+
+ |
| | | Added Matlab interface, which allows simulation
+of a decision rule in Matlab.
+
+ |
| | | Got rid of Matrix Template Library.
+
+ |
| | | Added checking of model residuals by the
+numerical integration. Three methods: checking along simulation path,
+checking along shocks, and on ellipse of states.
+
+ |
| | | Corrected a bug in calculation of higher moments
+of Normal dist.
+
+ |
| | | Corrected a bug of wrong drawing from Normal dist
+with non-zero covariances.
+
+ |
| |
+ | Added numerical integration module. Product and Smolyak
+quadratures over Gauss-Hermite and Gauss-Legendre, and quasi Monte
+Carlo.
+
+ |
+| 152
+ | 1.1
+ | 2005/04/22
+ |
+
+ |
| |
+ | Added a calculation of approximation at a stochastic steady state
+(still experimental).
+
+ |
| |
+ | Corrected a bug in Cholesky decomposition of variance-covariance
+matrix with off-diagonal elements.
+
+ |
+| 89
+ | 1.01
+ | 2005/02/23
+ |
+
+ |
| |
+ | Added version printout.
+
+ |
| |
+ | Corrected the bug of multithreading support for P4 HT processors running on Win32.
+
+ |
| |
+ | Enhanced Kronecker product code resulting in approx. 20% speedup.
+
+ |
| |
+ | Implemented vertical stack container refinement, and another
+method for sparse folded Faa Di Bruno (both not used yet).
+
+ |
+| 5
+ | 1.0
+ | 2005/02/23
+ | The first released version.
+
+ |
+
+
diff --git a/dynare++/doc/compiling-notes.txt b/dynare++/doc/compiling-notes.txt
new file mode 100644
index 000000000..275b76d10
--- /dev/null
+++ b/dynare++/doc/compiling-notes.txt
@@ -0,0 +1,23 @@
+It is suggested that you compile Dynare++ with gcc version 3.4. If you have
+other versions of gcc installed on your system, you need to select version
+3.4, for example,
+
+$ make "CC=gcc-3.4"
+
+For linking, you need to compile the required linear algebra libraries
+(blas, atlas, etc). Alternatively, you can install precompiled versions of
+these -- make sure that you select the version that matches your CPU.
+
+For example, if you have an SSE2 capable CPU, then on Debian GNU/Linux
+(etch) you need to install the following packages to get the precompiled
+linear algebra libraries:
+
+lapack3-dev
+atlas3-sse2-dev
+
+Then set the environment variable LD_LIBRARY_PATH, eg
+
+$ export LD_LIBRARY_PATH /usr/lib/sse2
+
+before calling make. This will include the shared libraries in the search
+path for ld.
diff --git a/dynare++/doc/dynare++-ramsey.tex b/dynare++/doc/dynare++-ramsey.tex
new file mode 100755
index 000000000..0044cb2d9
--- /dev/null
+++ b/dynare++/doc/dynare++-ramsey.tex
@@ -0,0 +1,157 @@
+\documentclass[10pt]{article}
+\usepackage{array,natbib,times}
+\usepackage{amsmath, amsthm, amssymb}
+
+%\usepackage[pdftex,colorlinks]{hyperref}
+
+\begin{document}
+
+\title{Implementation of Ramsey Optimal Policy in Dynare++, Timeless Perspective}
+
+\author{Ondra Kamen\'\i k}
+
+\date{June 2006}
+\maketitle
+
+\textbf{Abstract:} This document provides a derivation of Ramsey
+optimal policy from timeless perspective and describes its
+implementation in Dynare++.
+
+\section{Derivation of the First Order Conditions}
+
+Let us start with an economy populated by agents who take a number of
+variables exogenously, or given. These may include taxes or interest
+rates for example. These variables can be understood as decision (or control)
+variables of the timeless Ramsey policy (or social planner). The agent's
+information set at time $t$ includes mass-point distributions of these
+variables for all times after $t$. If $i_t$ denotes an interest rate
+for example, then the information set $I_t$ includes
+$i_{t|t},i_{t+1|t},\ldots,i_{t+k|t},\ldots$ as numbers. In addition
+the information set includes all realizations of past exogenous
+innovations $u_\tau$ for $\tau=t,t-1,\ldots$ and distibutions
+$u_\tau\sim N(0,\Sigma)$ for $\tau=t+1,\ldots$. These information sets will be denoted $I_t$.
+
+An information set including only the information on past realizations
+of $u_\tau$ and future distributions of $u_\tau\sim N(0\sigma)$ will
+be denoted $J_t$. We will use the following notation for expectations
+through these sets:
+\begin{eqnarray*}
+E^I_t[X] &=& E(X|I_t)\\
+E^J_t[X] &=& E(X|J_t)
+\end{eqnarray*}
+
+The agents optimize taking the decision variables of the social
+planner at $t$ and future as given. This means that all expectations
+they form are conditioned on the set $I_t$. Let $y_t$ denote a vector
+of all endogenous variables including the planer's decision
+variables. Let the number of endogenous variables be $n$. The economy
+can be described by $m$ equations including the first order conditions
+and transition equations:
+\begin{equation}\label{constr}
+E_t^I\left[f(y_{t-1},y_t,y_{t+1},u_t)\right] = 0.
+\end{equation}
+This lefts $n-m$
+the planner's control variables. The solution of this problem is a
+decision rule of the form:
+\begin{equation}\label{agent_dr}
+y_t=g(y_{t-1},u_t,c_{t|t},c_{t+1|t},\ldots,c_{t+k|t},\ldots),
+\end{equation}
+where $c$ is a vector of planner's control variables.
+
+Each period the social planner chooses the vector $c_t$ to maximize
+his objective such that \eqref{agent_dr} holds for all times following
+$t$. This would lead to $n-m$ first order conditions with respect to
+$c_t$. These first order conditions would contain unknown derivatives
+of endogenous variables with respect to $c$, which would have to be
+retrieved from the implicit constraints \eqref{constr} since the
+explicit form \eqref{agent_dr} is not known.
+
+The other way to proceed is to assume that the planner is so dumb that
+he is not sure what are his control variables. So he optimizes with
+respect to all $y_t$ given the constraints \eqref{constr}. If the
+planner's objective is $b(y_{t-1},y_t,y_{t+1},u_t)$ with a discount rate
+$\beta$, then the optimization problem looks as follows:
+\begin{align}
+\max_{\left\{y_\tau\right\}^\infty_t}&E_t^J
+\left[\sum_{\tau=t}^\infty\beta^{\tau-t}b(y_{\tau-1},y_\tau,y_{\tau+1},u_\tau)\right]\notag\\
+&\rm{s.t.}\label{planner_optim}\\
+&\hskip1cm E^I_\tau\left[f(y_{\tau-1},y_\tau,y_{\tau+1},u_\tau)\right]=0\quad\rm{for\ }
+\tau=\ldots,t-1,t,t+1,\ldots\notag
+\end{align}
+Note two things: First, each constraint \eqref{constr} in
+\eqref{planner_optim} is conditioned on $I_\tau$ not $I_t$. This is
+very important, since the behaviour of agents at period $\tau=t+k$ is
+governed by the constraint using expectations conditioned on $t+k$,
+not $t$. The social planner knows that at $t+k$ the agents will use
+all information available at $t+k$. Second, the constraints for the
+planner's decision made at $t$ include also constraints for agent's
+behaviour prior to $t$. This is because the agent's decision rules are
+given in the implicit form \eqref{constr} and not in the explicit form
+\eqref{agent_dr}.
+
+Using Lagrange multipliers, this can be rewritten as
+\begin{align}
+\max_{y_t}E_t^J&\left[\sum_{\tau=t}^\infty\beta^{\tau-t}b(y_{\tau-1},y_\tau,y_{\tau+1},u_\tau)\right.\notag\\
+&\left.+\sum_{\tau=-\infty}^{\infty}\beta^{\tau-t}\lambda^T_\tau E_\tau^I\left[f(y_{\tau-1},y_\tau,y_{\tau+1},u_\tau)\right]\right],
+\label{planner_optim_l}
+\end{align}
+where $\lambda_t$ is a vector of Lagrange multipliers corresponding to
+constraints \eqref{constr}. Note that the multipliers are multiplied
+by powers of $\beta$ in order to make them stationary. Taking a
+derivative wrt $y_t$ and putting it to zero yields the first order
+conditions of the planner's problem:
+\begin{align}
+E^J_t\left[\vphantom{\frac{\int^(_)}{\int^(\_)}}\right.&\frac{\partial}{\partial y_t}b(y_{t-1},y_t,y_{t+1},u_t)+
+\beta L^{+1}\frac{\partial}{\partial y_{t-1}}b(y_{t-1},y_t,y_{t+1},u_t)\notag\\
+&+\beta^{-1}\lambda_{t-1}^TE^I_{t-1}\left[L^{-1}\frac{\partial}{\partial y_{t+1}}f(y_{t-1},y_t,y_{t+1},u_t)\right]\notag\\
+&+\lambda_t^TE^I_t\left[\frac{\partial}{\partial y_{t}}f(y_{t-1},y_t,y_{t+1},u_t)\right]\notag\\
+&+\beta\lambda_{t+1}^TE^I_{t+1}\left[L^{+1}\frac{\partial}{\partial y_{t-1}}f(y_{t-1},y_t,y_{t+1},u_t)\right]
+\left.\vphantom{\frac{\int^(_)}{\int^(\_)}}\right]
+ = 0,\label{planner_optim_foc}
+\end{align}
+where $L^{+1}$ and $L^{-1}$ are one period lead and lag operators respectively.
+
+Now we have to make a few assertions concerning expectations
+conditioned on the different information sets to simplify
+\eqref{planner_optim_foc}. Recall the formula for integration through
+information on which another expectation is conditioned, this is:
+$$E\left[E\left[u|v\right]\right] = E[u],$$
+where the outer expectation integrates through $v$. Since $J_t\subset
+I_t$, by easy application of the above formula we obtain
+\begin{eqnarray}
+E^J_t\left[E^I_t\left[X\right]\right] &=& E^J_t\left[X\right]\quad\rm{and}\notag\\
+E^J_t\left[E^I_{t-1}\left[X\right]\right] &=& E^J_t\left[X\right]\label{e_iden}\\
+E^J_t\left[E^I_{t+1}\left[X\right]\right] &=& E^J_{t+1}\left[X\right]\notag
+\end{eqnarray}
+Now, the last term of \eqref{planner_optim_foc} needs a special
+attention. It is equal to
+$E^J_t\left[\beta\lambda^T_{t+1}E^I_{t+1}[X]\right]$. If we assume
+that the problem \eqref{planner_optim} has a solution, then there is a
+deterministic function from $J_{t+1}$ to $\lambda_{t+1}$ and so
+$\lambda_{t+1}\in J_{t+1}\subset I_{t+1}$. And the last term is equal
+to $E^J_{t}\left[E^I_{t+1}[\beta\lambda^T_{t+1}X]\right]$, which is
+$E^J_{t+1}\left[\beta\lambda^T_{t+1}X\right]$. This term can be
+equivalently written as
+$E^J_{t}\left[\beta\lambda^T_{t+1}E^J_{t+1}[X]\right]$. The reason why
+we write the term in this way will be clear later. All in all, we have
+\begin{align}
+E^J_t\left[\vphantom{\frac{\int^(_)}{\int^(\_)}}\right.&\frac{\partial}{\partial y_t}b(y_{t-1},y_t,y_{t+1},u_t)+
+\beta L^{+1}\frac{\partial}{\partial y_{t-1}}b(y_{t-1},y_t,y_{t+1},u_t)\notag\\
+&+\beta^{-1}\lambda_{t-1}^TL^{-1}\frac{\partial}{\partial y_{t+1}}f(y_{t-1},y_t,y_{t+1},u_t)\notag\\
+&+\lambda_t^T\frac{\partial}{\partial y_{t}}f(y_{t-1},y_t,y_{t+1},u_t)\notag\\
+&+\beta\lambda_{t+1}^TE^J_{t+1}\left[L^{+1}\frac{\partial}{\partial y_{t-1}}f(y_{t-1},y_t,y_{t+1},u_t)\right]
+\left.\vphantom{\frac{\int^(_)}{\int^(\_)}}\right]
+ = 0.\label{planner_optim_foc2}
+\end{align}
+Note that we have not proved that \eqref{planner_optim_foc} and
+\eqref{planner_optim_foc2} are equivalent. We proved only that if
+\eqref{planner_optim_foc} has a solution, then
+\eqref{planner_optim_foc2} is equivalent (and has the same solution).
+
+\section{Implementation}
+
+The user inputs $b(y_{t-1},y_t,y_{t+1},u_t)$, $\beta$, and agent's
+first order conditions \eqref{constr}. The algorithm has to produce
+\eqref{planner_optim_foc2}.
+
+\end{document}
diff --git a/dynare++/doc/dynare++-tutorial.tex b/dynare++/doc/dynare++-tutorial.tex
new file mode 100644
index 000000000..979ab50e1
--- /dev/null
+++ b/dynare++/doc/dynare++-tutorial.tex
@@ -0,0 +1,1512 @@
+\documentclass[10pt]{article}
+\usepackage{array,natbib}
+\usepackage{amsmath, amsthm, amssymb}
+
+\usepackage[pdftex,colorlinks]{hyperref}
+
+\begin{document}
+
+\title{DSGE Models with Dynare++. A Tutorial.}
+
+\author{Ondra Kamen\'\i k}
+
+\date{Jan 2009 (describes v. 1.3.7)}
+\maketitle
+
+\tableofcontents
+
+\section{Installation}
+
+Dynare++ installation procedure is pretty straightforward. Take the
+following steps:
+\begin{enumerate}
+\item Grab Dynare++ package from Dynare site either for Windows or for
+ Linux according to your operating system.\footnote{If unsure,
+ download the one for Windows.}
+\item Unzip or untar the package to a directory of your choice.
+\item Set operating system path to point to {\tt
+ dynare++-*.*}\footnote{Stars stand for a version number.}
+ subdirectory of the directory you have chosen. In Windows, this step
+ ensures that libraries distributed along Dynare++ are on the
+ path. This step is not really necessary in Linux.
+\item If you have Matlab and want to run custom simulations (see
+ \ref{custom}), then set the Matlab path to the path from the
+ previous step.
+\end{enumerate}
+
+Dynare++ uninstall procedure is even more simple. Just remove the directory {\tt
+ dynare++-*.*}.
+
+If you want (or need) to compile Dynare++ from sources, grab Dynare++
+source package and do your best.\footnote{Feel free to contact me if
+setting up the sources takes more time than solving your model by
+hand.}
+
+\section{Sample Session}
+
+As an example, let us take a simple DSGE model whose dynamic
+equilibrium is described by the following first order conditions:
+
+\begin{align*}
+&c_t\theta h_t^{1+\psi} = (1-\alpha)y_t\cr
+&\beta E_t\left[\frac{\exp(b_t)c_t}{\exp(b_{t+1})c_{t+1}}
+\left(\exp(b_{t+1})\alpha\frac{y_{t+1}}{k_{t+1}}+1-\delta\right)\right]=1\cr
+&y_t=\exp(a_t)k_t^\alpha h_t^{1-\alpha}\cr
+&k_{t}=\exp(b_{t-1})(y_{t-1}-c_{t-1})+(1-\delta)k_{t-1}\cr
+&a_t=\rho a_{t-1}+\tau b_{t-1}+\epsilon_t\cr
+&b_t=\tau a_{t-1}+\rho b_{t-1}+\nu_t
+\end{align*}
+
+\label{timing}
+The timing of this model is that the exogenous shocks $\epsilon_t$,
+and $\nu_t$ are observed by agents in the beginning of period $t$ and
+before the end of period $t$ all endogenous variables with index $t$
+are decided. The expectation operator $E_t$ works over the information
+accumulated just before the end of the period $t$ (this includes
+$\epsilon_t$, $\nu_t$ and all endogenous variables with index $t$).
+
+The exogenous shocks $\epsilon_t$ and $\nu_t$ are supposed to be
+serially uncorrelated with zero means and time-invariant
+variance-covariance matrix. In Dynare++, these variables are called
+exogenous; all other variables are endogenous. Now we are prepared to
+start writing a model file for Dynare++, which is an ordinary text
+file and could be created with any text editor.
+
+The model file starts with a preamble declaring endogenous and
+exogenous variables, parameters, and setting values of the
+parameters. Note that one can put expression on right hand sides. The
+preamble follows:
+
+{\small
+\begin{verbatim}
+var Y, C, K, A, H, B;
+varexo EPS, NU;
+
+parameters beta, rho, beta, alpha, delta, theta, psi, tau;
+alpha = 0.36;
+rho = 0.95;
+tau = 0.025;
+beta = 1/(1.03^0.25);
+delta = 0.025;
+psi = 0;
+theta = 2.95;
+\end{verbatim}
+}
+
+The section setting values of the parameters is terminated by a
+beginning of the {\tt model} section, which states all the dynamic
+equations. A timing convention of a Dynare++ model is the same as the
+timing of our example model, so we may proceed with writing the model
+equations. The time indexes of $c_{t-1}$, $c_t$, and $c_{t+1}$ are
+written as {\tt C(-1)}, {\tt C}, and {\tt C(1)} resp. The {\tt model}
+section looks as follows:
+
+{\small
+\begin{verbatim}
+model;
+C*theta*H^(1+psi) = (1-alpha)*Y;
+beta*exp(B)*C/exp(B(1))/C(1)*
+ (exp(B(1))*alpha*Y(1)/K(1)+1-delta) = 1;
+Y = exp(A)*K^alpha*H^(1-alpha);
+K = exp(B(-1))*(Y(-1)-C(-1)) + (1-delta)*K(-1);
+A = rho*A(-1) + tau*B(-1) + EPS;
+B = tau*A(-1) + rho*B(-1) + NU;
+end;
+\end{verbatim}
+}
+
+At this point, almost all information that Dynare++ needs has been
+provided. Only three things remain to be specified: initial values of
+endogenous variables for non-linear solver, variance-covariance matrix
+of the exogenous shocks and order of the Taylor approximation. Since
+the model is very simple, there is a closed form solution for the
+deterministic steady state. We use it as initial values for the
+non-linear solver. Note that the expressions on the right hand-sides in
+{\tt initval} section can reference values previously calculated. The
+remaining portion of the model file looks as follows:
+
+{\small
+\begin{verbatim}
+initval;
+A = 0;
+B = 0;
+H = ((1-alpha)/(theta*(1-(delta*alpha)
+ /(1/beta-1+delta))))^(1/(1+psi));
+Y = (alpha/(1/beta-1+delta))^(alpha/(1-alpha))*H;
+K = alpha/(1/beta-1+delta)*Y;
+C = Y - delta*K;
+end;
+
+vcov = [
+ 0.0002 0.00005;
+ 0.00005 0.0001
+];
+
+order = 7;
+\end{verbatim}
+}
+
+Note that the order of rows/columns of the variance-covariance matrix
+corresponds to the ordering of exogenous variables in the {\tt varexo}
+declaration. Since the {\tt EPS} was declared first, its variance is
+$0.0002$, and the variance of {\tt NU} is $0.0001$.
+
+Let the model file be saved as {\tt example1.mod}. Now we are prepared
+to solve the model. At the operating system command
+prompt\footnote{Under Windows it is a {\tt cmd} program, under Unix it
+is any shell} we issue a command:
+
+{\small
+\begin{verbatim}
+dynare++ example1.mod
+\end{verbatim}
+}
+
+When the program is finished, it produces two output files: a journal
+file {\tt example1.jnl} and a Matlab MAT-4 {\tt example1.mat}. The
+journal file contains information about time, memory and processor
+resources needed for all steps of solution. The output file is more
+interesting. It contains various simulation results. It can be loaded
+into Matlab or Scilab and examined.%
+\footnote{For Matlab {\tt load example1.mat}, for Scilab {\tt
+mtlb\_load example1.mat}} The following examples are done in Matlab,
+everything would be very similar in Scilab.
+
+Let us first examine the contents of the MAT file:
+{\small
+\begin{verbatim}
+>> load example1.mat
+>> who
+
+Your variables are:
+
+dyn_g_1 dyn_i_Y dyn_npred
+dyn_g_2 dyn_irfm_EPS_mean dyn_nstat
+dyn_g_3 dyn_irfm_EPS_var dyn_shocks
+dyn_g_4 dyn_irfm_NU_mean dyn_ss
+dyn_g_5 dyn_irfm_NU_var dyn_state_vars
+dyn_i_A dyn_irfp_EPS_mean dyn_steady_states
+dyn_i_B dyn_irfp_EPS_var dyn_vars
+dyn_i_C dyn_irfp_NU_mean dyn_vcov
+dyn_i_EPS dyn_irfp_NU_var dyn_vcov_exo
+dyn_i_H dyn_mean
+dyn_i_K dyn_nboth
+dyn_i_NU dyn_nforw
+\end{verbatim}
+}
+
+All the variables coming from one MAT file have a common prefix. In
+this case it is {\tt dyn}, which is Dynare++ default. The prefix can
+be changed, so that the multiple results could be loaded into one Matlab
+session.
+
+In the default setup, Dynare++ solves the Taylor approximation to the
+decision rule and calculates unconditional mean and covariance of the
+endogenous variables, and generates impulse response functions. The
+mean and covariance are stored in {\tt dyn\_mean} and {\tt
+dyn\_vcov}. The ordering of the endogenous variables is given by {\tt
+dyn\_vars}.
+
+In our example, the ordering is
+
+{\small
+\begin{verbatim}
+>> dyn_vars
+dyn_vars =
+H
+A
+Y
+C
+K
+B
+\end{verbatim}
+}
+
+and unconditional mean and covariance are
+
+{\small
+\begin{verbatim}
+>> dyn_mean
+dyn_mean =
+ 0.2924
+ 0.0019
+ 1.0930
+ 0.8095
+ 11.2549
+ 0.0011
+>> dyn_vcov
+dyn_vcov =
+ 0.0003 0.0006 0.0016 0.0004 0.0060 0.0004
+ 0.0006 0.0024 0.0059 0.0026 0.0504 0.0012
+ 0.0016 0.0059 0.0155 0.0069 0.1438 0.0037
+ 0.0004 0.0026 0.0069 0.0040 0.0896 0.0016
+ 0.0060 0.0504 0.1438 0.0896 2.1209 0.0405
+ 0.0004 0.0012 0.0037 0.0016 0.0405 0.0014
+\end{verbatim}
+}
+
+The ordering of the variables is also given by indexes starting with
+{\tt dyn\_i\_}. Thus the mean of capital can be retrieved as
+
+{\small
+\begin{verbatim}
+>> dyn_mean(dyn_i_K)
+ans =
+ 11.2549
+\end{verbatim}
+}
+
+\noindent and covariance of labor and capital by
+
+{\small
+\begin{verbatim}
+>> dyn_vcov(dyn_i_K,dyn_i_H)
+ans =
+ 0.0060
+\end{verbatim}
+}
+
+The impulse response functions are stored in matrices as follows
+\begin{center}
+\begin{tabular}{|l|l|}
+\hline
+matrix& response to\\
+\hline
+{\tt dyn\_irfp\_EPS\_mean}& positive impulse to {\tt EPS}\\
+{\tt dyn\_irfm\_EPS\_mean}& negative impulse to {\tt EPS}\\
+{\tt dyn\_irfp\_NU\_mean}& positive impulse to {\tt NU}\\
+{\tt dyn\_irfm\_NU\_mean}& negative impulse to {\tt NU}\\
+\hline
+\end{tabular}
+\end{center}
+All shocks sizes are one standard error. Rows of the matrices
+correspond to endogenous variables, columns correspond to
+periods. Thus capital response to a positive shock to {\tt EPS} can be
+plotted as
+
+{\small
+\begin{verbatim}
+plot(dyn_irfp_EPS_mean(dyn_i_K,:));
+\end{verbatim}
+}
+
+The data is in units of the respective variables, so in order to plot
+the capital response in percentage changes from the decision rule's
+fix point (which is a vector {\tt dyn\_ss}), one has to issue the
+commands:
+
+{\small
+\begin{verbatim}
+Kss=dyn_ss(dyn_i_K);
+plot(100*dyn_irfp_EPS_mean(dyn_i_K,:)/Kss);
+\end{verbatim}
+}
+
+The plotted impulse response shows that the model is pretty persistent
+and that the Dynare++ default for a number of simulated periods is not
+sufficient. In addition, the model persistence puts in doubt also a
+number of simulations. The Dynare++ defaults can be changed when
+calling Dynare++, in operating system's command prompt, we issue a
+command:
+
+{\small
+\begin{verbatim}
+dynare++ --per 300 --sim 150 example1.mod
+\end{verbatim}
+}
+
+\noindent This sets the number of simulations to $150$ and the number
+of periods to $300$ for each simulation giving $45000$ total simulated
+periods.
+
+\section{Sample Optimal Policy Session}
+\label{optim_tut}
+
+Suppose that one wants to solve the following optimal policy problem
+with timeless perspective.\footnote{See \ref{ramsey} on how to solve
+Ramsey optimality problem within this framework} The following
+optimization problem is how to choose capital taxes financing public
+good to maximize agent's utility from consumption good and public
+good. The problem takes the form:
+\begin{align*}
+\max_{\{\tau_t\}_{t_0}^\infty}
+E_{t_0}\sum_{t=t_0}^\infty &\beta^{t-t_0}\left(u(c_t)+av(g_t)\right)\\
+\hbox{subject\ to}&\\
+u'(c_t) &=
+\beta E_t\left[u'(c_{t+1})\left(1-\delta+f'(k_{t+1})(1-\alpha\tau_{t+1})\right)\right]\\
+K_t &= (1-\delta)K_{t-1} + (f(K_{t-1}) - c_{t-1} - g_{t-1})\\
+g_t &= \tau_t\alpha f(K_t),\\
+\hbox{where\ } t & = \ldots,t_0-1,t_0,t_0+1,\ldots
+\end{align*}
+$u(c_t)$ is utility from consuming the consumption good, $v(g_t)$ is
+utility from consuming the public good, $f(K_t)$ is a production
+function $f(K_t) = Z_tK_t^\alpha$. $Z_t$ is a technology shock modeled
+as AR(1) process. The three constraints come from the first order
+conditions of a representative agent. We suppose that it pursues a
+different objective, namely lifetime utility involving only
+consumption $c_t$. The representative agents chooses between
+consumption and investment. It rents the capital to firms and supplies
+constant amount of labour. All output is paid back to consumer in form
+of wage and capital rent. Only the latter is taxed. We suppose that
+the optimal choice has been taking place from infinite past and will
+be taking place for ever. Further we suppose the same about the
+constraints.
+
+Let us choose the following functional forms:
+\begin{eqnarray*}
+u(c_t) &=& \frac{c_t^{1-\eta}}{1-\eta}\\
+v(g_t) &=& \frac{g_t^{1-\phi}}{1-\phi}\\
+f(K_t) &=& K_t^\alpha
+\end{eqnarray*}
+
+Then the problem can be coded into Dynare++ as follows. We start with
+a preamble which states all the variables, shocks and parameters:
+{\small
+\begin{verbatim}
+var C G K TAU Z;
+
+varexo EPS;
+
+parameters eta beta alpha delta phi a rho;
+
+eta = 2;
+beta = 0.99;
+alpha = 0.3;
+delta = 0.10;
+phi = 2.5;
+a = 0.1;
+rho = 0.7;
+\end{verbatim}
+}
+
+Then we specify the planner's objective and the discount factor in the
+objective. The objective is an expression (possibly including also
+variable leads and lags), and the discount factor must be one single
+declared parameter:
+{\small
+\begin{verbatim}
+planner_objective C^(1-eta)/(1-eta) + a*G^(1-phi)/(1-phi);
+
+planner_discount beta;
+\end{verbatim}
+}
+
+The model section will contain only the constraints of the social
+planner. These are capital accumulation, identity for the public
+product, AR(1) process for $Z_t$ and the first order condition of the
+representative agent (with different objective).
+{\small
+\begin{verbatim}
+model;
+K = (1-delta)*K(-1) + (exp(Z(-1))*K(-1)^alpha - C(-1) - G(-1));
+G = TAU*alpha*K^alpha;
+Z = rho*Z(-1) + EPS;
+C^(-eta) = beta*C(+1)^(-eta)*(1-delta +
+ exp(Z(+1))*alpha*K(+1)^(alpha-1)*(1-alpha*TAU(+1)));
+end;
+\end{verbatim}
+}
+
+Now we have to provide a good guess for non-linear solver calculating
+the deterministic steady state. The model's steady state has a closed
+form solution if the taxes are known. So we provide a guess for
+taxation {\tt TAU} and then use the closed form solution for capital,
+public good and consumption:\footnote{Initial guess for Lagrange
+multipliers and some auxiliary variables is calculated automatically. See
+\ref{opt_init} for more details.}
+{\small
+\begin{verbatim}
+initval;
+TAU = 0.70;
+K = ((delta+1/beta-1)/(alpha*(1-alpha*TAU)))^(1/(alpha-1));
+G = TAU*alpha*K^alpha;
+C = K^alpha - delta*K - G;
+Z = 0;
+\end{verbatim}
+}
+
+Finally, we have to provide the order of approximation, and the
+variance-covariance matrix of the shocks (in our case we have only one
+shock):
+{\small
+\begin{verbatim}
+order = 4;
+
+vcov = [
+ 0.01
+];
+\end{verbatim}
+}
+
+After this model file has been run, we can load the resulting MAT-file
+into the Matlab (or Scilab) and examine its contents:
+{\small
+\begin{verbatim}
+>> load kp1980_2.mat
+>> who
+
+Your variables are:
+
+dyn_g_1 dyn_i_MULT1 dyn_nforw
+dyn_g_2 dyn_i_MULT2 dyn_npred
+dyn_g_3 dyn_i_MULT3 dyn_nstat
+dyn_g_4 dyn_i_TAU dyn_shocks
+dyn_i_AUX_3_0_1 dyn_i_Z dyn_ss
+dyn_i_AUX_4_0_1 dyn_irfm_EPS_mean dyn_state_vars
+dyn_i_C dyn_irfm_EPS_var dyn_steady_states
+dyn_i_EPS dyn_irfp_EPS_mean dyn_vars
+dyn_i_G dyn_irfp_EPS_var dyn_vcov
+dyn_i_K dyn_mean dyn_vcov_exo
+dyn_i_MULT0 dyn_nboth
+\end{verbatim}
+}
+
+The data dumped into the MAT-file have the same structure as in the
+previous example of this tutorial. The only difference is that
+Dynare++ added a few more variables. Indeed:
+{\small
+\begin{verbatim}
+>> dyn_vars
+dyn_vars =
+MULT1
+G
+MULT3
+C
+K
+Z
+TAU
+AUX_3_0_1
+AUX_4_0_1
+MULT0
+MULT2
+\end{verbatim}
+}
+Besides the five variables declared in the model ({\tt C}, {\tt G},
+{\tt K}, {\tt TAU}, and {\tt Z}), Dy\-na\-re++ added 6 more, four as Lagrange
+multipliers of the four constraints, two as auxiliary variables for
+shifting in time. See \ref{aux_var} for more details.
+
+The structure and the logic of the MAT-file is the same as these new 6
+variables were declared in the model file and the file is examined in
+the same way.
+
+For instance, let us examine the Lagrange multiplier of the optimal
+policy associated with the consumption first order condition. Recall
+that the consumers' objective is different from the policy
+objective. Therefore, the constraint will be binding and the
+multiplier will be non-zero. Indeed, its deterministic steady state,
+fix point and mean are as follows:
+{\small
+\begin{verbatim}
+>> dyn_steady_states(dyn_i_MULT3,1)
+ans =
+ -1.3400
+>> dyn_ss(dyn_i_MULT3)
+ans =
+ -1.3035
+>> dyn_mean(dyn_i_MULT3)
+ans =
+ -1.3422
+\end{verbatim}
+}
+
+\section{What Dynare++ Calculates}
+\label{dynpp_calc}
+
+Dynare++ solves first order conditions of a DSGE model in the recursive form:
+\begin{equation}\label{focs}
+E_t[f(y^{**}_{t+1},y_t,y^*_{t-1},u_t)]=0,
+\end{equation}
+where $y$ is a vector of endogenous variables, and $u$ a vector of
+exogenous variables. Some of elements of $y$ can occur at time $t+1$,
+these are $y^{**}$. Elements of $y$ occurring at time $t-1$ are denoted
+$y^*$. The exogenous shocks are supposed to be serially independent
+and normally distributed $u_t\sim N(0,\Sigma)$.
+
+The solution of this dynamic system is a decision rule
+\[
+y_t=g(y^*_{t-1},u_t)
+\]
+Dynare++ calculates a Taylor approximation of this decision rule of a
+given order. The approximation takes into account deterministic
+effects of future volatility, so a point about which the Taylor
+approximation is done will be different from the fix point $y$ of the rule
+yielding $y=g(y^*,0)$.
+
+The fix point of a rule corresponding to a model with $\Sigma=0$ is
+called {\it deterministic steady state} denoted as $\bar y$. In
+contrast to deterministic steady state, there is no consensus in
+literature how to call a fix point of the rule corresponding to a
+model with non-zero $\Sigma$. I am tempted to call it {\it stochastic
+ steady state}, however, it might be confused with unconditional mean
+or with steady distribution. So I will use a term {\it fix point} to
+avoid a confusion.
+
+By default, Dynare++ solves the Taylor approximation about the
+deterministic steady state. Alternatively, Dynare++ can split the
+uncertainty to a few steps and take smaller steps when calculating the
+fix points. This is controlled by an option {\tt --steps}. For the
+brief description of the second method, see \ref{multistep_alg}.
+
+\subsection{Decision Rule Form}
+\label{dr_form}
+
+In case of default solution algorithm (approximation about the
+deterministic steady state $\bar y$), Dynare++ calculates the higher
+order derivatives of the equilibrium rule to get a decision rule of
+the following form. In Einstein notation, it is:
+\[
+y_t-\bar y = \sum_{i=0}^k\frac{1}{i!}\left[g_{(y^*u)^i}\right]
+_{\alpha_1\ldots\alpha_i}
+\prod_{j=1}^i\left[\begin{array}{c} y^*_{t-1}-\bar y^*\\ u_t \end{array}\right]
+^{\alpha_j}
+\]
+
+Note that the ergodic mean will be different from the deterministic
+steady state $\bar y$ and thus deviations $y^*_{t-1}-\bar y^*$ will
+not be zero in average. This implies that in average we will commit
+larger round off errors than if we used the decision rule expressed in
+deviations from a point closer to the ergodic mean. Therefore, by
+default, Dynare++ recalculates this rule and expresses it in
+deviations from the stochastic fix point $y$.
+\[
+y_t-y = \sum_{i=1}^k\frac{1}{i!}\left[\tilde g_{(y^*u)^i}\right]
+_{\alpha_1\ldots\alpha_i}
+\prod_{j=1}^i\left[\begin{array}{c} y^*_{t-1}-y^*\\ u_t \end{array}\right]
+^{\alpha_j}
+\]
+Note that since the rule is centralized around its fix point, the
+first term (for $i=0$) drops out.
+
+Also note, that this rule mathematically equivalent to the rule
+expressed in deviations from the deterministic steady state, and still
+it is an approximation about the deterministic steady state. The fact
+that it is expressed in deviations from a different point should not
+be confused with the algorithm in \ref{multistep_alg}.
+
+This centralization can be avoided by invoking {\tt --no-centralize}
+command line option.
+
+\subsection{Taking Steps in Volatility Dimension}
+\label{multistep_alg}
+
+For models, where volatility of the exogenous shocks plays a big
+role, the approximation about deterministic steady state can be poor,
+since the equilibrium dynamics can be very different from the dynamics
+in the vicinity of the perfect foresight (deterministic steady state).
+
+Therefore, Dynare++ has on option {\tt --steps} triggering a multistep
+algorithm. The algorithm splits the volatility to a given number of
+steps. Dynare++ attempts to calculate approximations about fix points
+corresponding to these levels of volatility. The problem is that if we
+want to calculate higher order approximations about fix points
+corresponding to volatilities different from zero (as in the case of
+deterministic steady state), then the derivatives of lower orders
+depend on derivatives of higher orders with respect to forward looking
+variables. The multistep algorithm in each step approximates the
+missing higher order derivatives with extrapolations based on the
+previous step.
+
+In this way, the approximation of the stochastic fix point and the
+derivatives about this fix point are obtained. It is difficult to a
+priori decide whether this algorithm yields a better decision
+rule. Nothing is guaranteed, and the resulted decision rule should be
+checked with a numerical integration. See \ref{checks}.
+
+\subsection{Simulating the Decision Rule}
+
+After some form of a decision rule is calculated, it is simulated to
+obtain draws from ergodic (unconditional) distribution of endogenous
+variables. The mean and the covariance are reported. There are two
+ways how to calculate the mean and the covariance. The first one is to
+store all simulated samples and calculate the sample mean and
+covariance. The second one is to calculate mean and the covariance in
+the real-time not storing the simulated sample. The latter case is
+described below (see \ref{rt_simul}).
+
+The stored simulated samples are then used for impulse response
+function calculations. For each shock, the realized shocks in these
+simulated samples (control simulations) are taken and an impulse is
+added and the new realization of shocks is simulated. Then the control
+simulation is subtracted from the simulation with the impulse. This is
+done for all control simulations and the results are averaged. As the
+result, we get an expectation of difference between paths with impulse
+and without impulse. In addition, the sample variances are
+reported. They might be useful for confidence interval calculations.
+
+For each shock, Dynare++ calculates IRF for two impulses, positive and
+negative. Size of an impulse is one standard error of a respective
+shock.
+
+The rest of this subsection is divided to three parts giving account
+on real-time simulations, conditional simulations, and on the way how
+random numbers are generated resp.
+
+\subsubsection{Simulations With Real-Time Statistics}
+\label{rt_simul}
+
+When one needs to simulate large samples to get a good estimate of
+unconditional mean, simulating the decision rule with statistics
+calculated in real-time comes handy. The main reason is that the
+storing of all simulated samples may not fit into the available
+memory.
+
+The real-time statistics proceed as follows: We model the ergodic
+distribution as having normal distribution $y\sim N(\mu,\Sigma)$. Further,
+the parameters $\mu$ and $\Sigma$ are modelled as:
+\begin{eqnarray*}
+ \Sigma &\sim& {\rm InvWishart}_\nu(\Lambda)\\
+ \mu|\Sigma &\sim& N(\bar\mu,\Sigma/\kappa) \\
+\end{eqnarray*}
+This model of $p(\mu,\Sigma)$ has an advantage of conjugacy, i.e. a
+prior distribution has the same form as posterior. This property is
+used in the calculation of real-time estimates of $\mu$ and $\Sigma$,
+since it suffices to maintain only the parameters of $p(\mu,\Sigma)$
+conditional observed draws so far. The parameters are: $\nu$,
+$\Lambda$, $\kappa$, and $\bar\mu$.
+
+The mean of $\mu,\Sigma|Y$, where $Y$ are all the draws (simulated
+periods) is reported.
+
+\subsubsection{Conditional Distributions}
+\label{cond_dist}
+
+Starting with version 1.3.6, Dynare++ calculates variable
+distributions $y_t$ conditional on $y_0=\bar y$, where $\bar y$ is the
+deterministic steady state. If triggered, Dynare++ simulates a given
+number of samples with a given number of periods all starting at
+the deterministic steady state. Then for each time $t$, mean
+$E[y_t|y_0=\bar y]$ and variances $E[(y_t-E[y_t|y_0=\bar
+y])(y_t-E[y_t|y_0=\bar y])^T|y_0=\bar y]$ are reported.
+
+\subsubsection{Random Numbers}
+\label{random_numbers}
+
+For generating of the pseudo random numbers, Dynare++ uses Mersenne
+twister by Makoto Matsumoto and Takuji Nishimura. Because of the
+parallel nature of Dynare++ simulations, each simulated sample gets
+its own instance of the twister. Each such instance is seeded before
+the simulations are started. This is to prevent additional randomness
+implied by the operating system's thread scheduler to interfere with
+the pseudo random numbers.
+
+For seeding the individual instances of the Mersenne twister assigned
+to each simulated sample the system (C library) random generator is
+used. These random generators do not have usually very good
+properties, but we use them only to seed the Mersenne twister
+instances. The user can set the initial seed of the system random
+generator and in this way deterministically choose the seeds of all
+instances of the Mersenne twister.
+
+In this way, it is guaranteed that two runs of Dynare++
+with the same seed will yield the same results regardless the
+operating system's scheduler. The only difference may be caused by a
+different round-off errors committed when the same set of samples are
+summed in the different order (due to the operating system's scheduler).
+
+\subsection{Numerical Approximation Checks}
+\label{checks}
+
+Optionally, Dynare++ can run three kinds of checks for Taylor
+approximation errors. All three methods numerically calculate
+the residual of the DSGE equations
+\[
+E[f(g^{**}(g^*(y^*,u),u'),g(y^*,u),y^*,u)|y^*,u]
+\]
+which must be ideally zero for all $y^*$ and $u$. This integral is
+evaluated by either product or Smolyak rule applied to one dimensional
+Gauss--Hermite quadrature. The user does not need to care about the
+decision. An algorithm yielding higher quadrature level and less
+number of evaluations less than a user given maximum is selected.
+
+The three methods differ only by a set of $y^*$ and $u$ where the
+residuals are evaluated. These are:
+\begin{itemize}
+\item The first method calculates the residuals along the shocks for
+fixed $y^*$ equal to the fix point. We let all elements of $u$ be
+fixed at $0$ but one element, which varies from $-\mu\sigma$ to
+$\mu\sigma$, where $\sigma$ is a standard error of the element and
+$\mu$ is the user given multiplier. In this way we can see how the
+approximation error grows if the fix point is disturbed by a shock of
+varying size.
+\item The second method calculates the residuals along a simulation
+path. A random simulation is run, and at each point the residuals are
+reported.
+\item The third method calculates the errors on an ellipse of the
+state variables $y^*$. The shocks $u$ are always zero. The ellipse is
+defined as
+\[\{Ax|\; \Vert x\Vert_2=\mu\},\]
+where $\mu$ is a user given multiplier, and $AA^T=V$ for $V$ being a
+covariance of endogenous variables based on the first order
+approximation. The method calculates the residuals at low discrepancy
+sequence of points on the ellipse. Both the residuals and the points
+are reported.
+\end{itemize}
+
+\section{Optimal Policy with Dynare++}
+\label{optim}
+
+Starting with version 1.3.2, Dynare++ is able to automatically
+generate and then solve the first order conditions for a given
+objective and (possibly) forward looking constraints. Since the
+constraints can be forward looking, the use of this feature will
+mainly be in optimal policy or control.
+
+The only extra thing which needs to be added to the model file is a
+specification of the policy's objective. This is done by two keywords,
+placed not before parameter settings. If the objective is to maximize
+$$E_{t_0}\sum_{t=t_0}^\infty\beta^{t-t_0}\left[\frac{c_t^{1-\eta}}{1-\eta}+
+a\frac{g_t^{1-\phi}}{1-\phi}\right],$$
+then the keywords will be:
+{\small
+\begin{verbatim}
+planner_objective C^(1-eta)/(1-eta) + a*G^(1-phi)/(1-phi);
+
+planner_discount beta;
+\end{verbatim}
+}
+
+Dynare++ parses the file and if the two keywords are present, it
+automatically derives the first order conditions for the problem. The
+first order conditions are put to the form \eqref{focs} and solved. In
+this case, the equations in the {\tt model} section are understood as
+the constraints (they might come as the first order conditions from
+optimizations of other agents) and their number must be less than the
+number of endogenous variables.
+
+This section further describes how the optimal policy first order
+conditions look like, then discusses some issues with the initial
+guess for deterministic steady state, and finally describes how to
+simulate Ramsey policy within this framework.
+
+\subsection{First Order Conditions}
+
+Mathematically, the optimization problem looks as follows:
+\begin{align}
+\max_{\left\{y_\tau\right\}^\infty_t}&E_t
+\left[\sum_{\tau=t}^\infty\beta^{\tau-t}b(y_{\tau-1},y_\tau,y_{\tau+1},u_\tau)\right]\notag\\
+&\rm{s.t.}\label{planner_optim}\\
+&\hskip1cm E^I_\tau\left[f(y_{\tau-1},y_\tau,y_{\tau+1},u_\tau)\right]=0\quad\rm{for\ }
+\tau=\ldots,t-1,t,t+1,\ldots\notag
+\end{align}
+where $E^I$ is an expectation operator over an information set including,
+besides all the past, all future realizations of policy's control
+variables and distributions of future shocks $u_t\sim
+N(0,\Sigma)$. The expectation operator $E$ integrates over an
+information including only distributions of $u_t$ (besides the past).
+
+Note that the constraints $f$ take place at all times, and they are
+conditioned at the running $\tau$ since the policy knows that the
+agents at time $\tau$ will use all the information available at
+$\tau$.
+
+The maximization problem can be rewritten using Lagrange multipliers as:
+\begin{align}
+\max_{y_t}E_t&\left[\sum_{\tau=t}^\infty\beta^{\tau-t}b(y_{\tau-1},y_\tau,y_{\tau+1},u_\tau)\right.\notag\\
+&\left.+\sum_{\tau=-\infty}^{\infty}\beta^{\tau-t}\lambda^T_\tau E_\tau^I\left[f(y_{\tau-1},y_\tau,y_{\tau+1},u_\tau)\right]\right],
+\label{planner_optim_l}
+\end{align}
+where $\lambda_t$ is a column vector of Lagrange multipliers.
+
+After some manipulations with compounded expectations over different
+information sets, one gets the following first order conditions:
+\begin{align}
+E_t\left[\vphantom{\frac{\int^(_)}{\int^(\_)}}\right.&\frac{\partial}{\partial y_t}b(y_{t-1},y_t,y_{t+1},u_t)+
+\beta L^{+1}\frac{\partial}{\partial y_{t-1}}b(y_{t-1},y_t,y_{t+1},u_t)\notag\\
+&+\beta^{-1}\lambda_{t-1}^TL^{-1}\frac{\partial}{\partial y_{t+1}}f(y_{t-1},y_t,y_{t+1},u_t)\notag\\
+&+\lambda_t^T\frac{\partial}{\partial y_{t}}f(y_{t-1},y_t,y_{t+1},u_t)\notag\\
+&+\beta\lambda_{t+1}^TE_{t+1}\left[L^{+1}\frac{\partial}{\partial y_{t-1}}f(y_{t-1},y_t,y_{t+1},u_t)\right]
+\left.\vphantom{\frac{\int^(_)}{\int^(\_)}}\right]
+ = 0,\label{planner_optim_foc2}
+\end{align}
+where $L^{+1}$ is one period lead operator, and $L^{-1}$ is one period lag operator.
+
+Dynare++ takes input corresponding to \eqref{planner_optim},
+introduces the Lagrange multipliers according to
+\eqref{planner_optim_l}, and using its symbolic derivator it compiles
+\eqref{planner_optim_foc2}. The system \eqref{planner_optim_foc2} with
+the constraints from \eqref{planner_optim_l} is then solved in the
+same way as the normal input \eqref{focs}.
+
+\subsection{Initial Guess for Deterministic Steady State}
+\label{opt_init}
+
+Solving deterministic steady state of non-linear dynamic systems is
+not trivial and the first order conditions for optimal policy add
+significant complexity. The {\tt initval} section allows to input the
+initial guess of the non-linear solver. It requires that all user
+declared endogenous variables be initialized. However, in most cases,
+we have no idea what are good initial guesses for the Lagrange
+multipliers.
+
+For this reason, Dynare++ calculates an initial guess of Lagrange
+multipliers using user provided initial guesses of all other
+endogenous variables. It uses the linearity of the Lagrange
+multipliers in the \eqref{planner_optim_foc2}. In its static form,
+\eqref{planner_optim_foc2} looks as follows:
+\begin{align}
+&\frac{\partial}{\partial y_t}b(y,y,y,0)+
+\beta\frac{\partial}{\partial y_{t-1}}b(y,y,y,0)\notag\\
+&+\lambda^T\left[\beta^{-1}\frac{\partial}{\partial y_{t+1}}f(y,y,y,0)
+ +\frac{\partial}{\partial y_{t}}f(y,y,y,0)
+ +\beta\frac{\partial}{\partial y_{t-1}}f(y,y,y,0)\right]
+ = 0\label{planner_optim_static}
+\end{align}
+
+The user is required to provide an initial guess of all declared
+variables (all $y$). Then \eqref{planner_optim_static} becomes an
+overdetermined linear system in $\lambda$, which is solved by means of
+the least squares. The closer the initial guess of $y$ is to the exact
+solution, the closer are the Lagrange multipliers $\lambda$.
+
+The calculated Lagrange multipliers by the least squares are not used,
+if they are set in the {\tt initval} section. In other words, if a
+multiplier has been given a value in the {\tt initval} section, then
+the value is used, otherwise the calculated value is taken.
+
+For even more difficult problems, Dynare++ generates two Matlab files
+calculating a residual of the static system and its derivative. These
+can be used in Matlab's {\tt fsolve} or other algorithm to get an
+exact solution of the deterministic steady state. See
+\ref{output_matlab_scripts} for more details.
+
+Finally, Dynare++ might generate a few auxiliary variables. These are
+simple transformations of other variables. They are initialized
+automatically and the user usually does not need to care about it.
+
+\subsection{Optimal Ramsey Policy}
+\label{ramsey}
+
+Dynare++ solves the optimal policy problem with timeless
+perspective. This means that it assumes that the constraints in
+\eqref{planner_optim} are valid from the infinite past to infinite
+future. Dynare++ calculation of ergodic distribution then assumes that
+the policy has been taking place from infinite past.
+
+If some constraints in \eqref{planner_optim} are forward looking, this
+will result in some backward looking Lagrange multipliers. Such
+multipliers imply possibly time inconsistent policy in the states of
+the ``original'' economy, since these backward looking multipliers add
+new states to the ``optimized'' economy. In this respect, the timeless
+perspective means that there is no fixed initial distribution of such
+multipliers, instead, their ergodic distribution is taken.
+
+In contrast, Ramsey optimal policy is started at $t=0$. This means
+that the first order conditions at $t=0$ are different than the first
+order conditions at $t\geq 1$, which are
+\eqref{planner_optim_foc2}. However, it is not difficult to assert
+that the first order conditions at $t=0$ are in the form of
+\eqref{planner_optim_foc2} if all the backward looking Lagrange
+multipliers are set to zeros at period $-1$, i.e. $\lambda_{-1}=0$.
+
+All in all, the solution of \eqref{planner_optim_foc2} calculated by
+Dynare++ can be used as a Ramsey optimal policy solution provided that
+all the backward looking Lagrange multipliers were set to zeros prior
+to the first simulation period. This can be done by setting the
+initial state of a simulation path in {\tt dynare\_simul.m}. If this
+is applied on the example from \ref{optim_tut}, then we may do the
+following in the command prompt:
+{\small
+\begin{verbatim}
+>> load kp1980_2.mat
+>> shocks = zeros(1,100);
+>> ystart = dyn_ss;
+>> ystart(dyn_i_MULT3) = 0;
+>> r=dynare_simul('kp1980_2.mat',shocks,ystart);
+\end{verbatim}
+}
+This will simulate the economy if the policy was introduced in the
+beginning and no shocks happened.
+
+More information on custom simulations can be obtained by typing:
+{\small
+\begin{verbatim}
+help dynare_simul
+\end{verbatim}
+}
+
+
+\section{Running Dynare++}
+
+This section deals with Dynare++ input. The first subsection
+\ref{dynpp_opts} provides a list of command line options, next
+subsection \ref{dynpp_mod} deals with a format of Dynare++ model file,
+and the last subsection discusses incompatibilities between Dynare
+Matlab and Dynare++.
+
+\subsection{Command Line Options}
+\label{dynpp_opts}
+
+The calling syntax of the Dynare++ is
+
+{\small
+\begin{verbatim}
+dynare++ [--help] [--version] [options]