\documentclass[11pt,a4paper]{article}
\usepackage{bibmad,graphicx,latexsym,amssymb,times}
\usepackage[cp850]{inputenc}

\begin{document}
\title{Stochastic simulations with {\sc Dynare}. \\ A practical guide.}
\author{Fabrice Collard (GREMAQ, University of Toulouse)\\Adapted for Dynare 4.1\\ by Michel Juillard and S\'ebastien Villemot (CEPREMAP)}
\date{First draft: February 2001\hspace{10mm}This draft: December 2009.}
\maketitle
This document describes a model involving both endogenous and exogenous state variable. We first describe the theoretical model, before showing how the perturbation method is implemented in {\sc Dynare}.

\section{A theoretical model} 
We consider an economy that consists of a large number of dynastic households and a large number of firms. Firms are producing a homogeneous final product that can be either consumed or invested by means of capital and labor services. Firms own their capital stock and hire labor supplied by the households. Households own the firms. In each and every period three perfectly competitive markets open --- the markets for consumption goods, labor services, and financial capital in the form of firms' shares. 
Household preferences are characterized by the lifetime utility function:
\begin{equation}
  E_t\sum_{\tau=t}^{\infty}{\beta^\star}^{\tau-t} \left(\log(c_t)-\theta\frac{h_t^{1+\psi}}{1+\psi}\right)
  \label{eq:ut}
\end{equation}
\noindent where $0<\beta^\star<1$ is a constant discount factor, $c_t$ is consumption in period
$t$, $h_t$ is the fraction of total available time devoted to productive activity in period $t$, $\theta>0$ and $\psi\geqslant 0$. We assume that there exists a central planner that determines hours, consumption and capital accumulation maximizing the household's utility function subject to the following budget constraint 
\begin{equation}
  c_t+i_t=y_t
  \label{eq:bud}
\end{equation}
\noindent where $i_t$ is investment and $y_t$ is output. Investment is used to form physical capital, which accumulates in the standard form as:
\begin{equation}
k_{t+1}=\exp(b_t) i_t+(1-\delta)k_t \mbox{ with } 0<\delta<1
\label{eq:acc}
\end{equation}
where $\delta$ is the constant physical depreciation rate. $b_t$ is a shock affecting incorporated technological progress, which properties will be defined later.

Output is produced by means of capital and labor services, relying on a constant returns to scale technology represented by the following Cobb--Douglas production function:
\begin{equation}
 y_t=\exp(a_t) k_t^\alpha h_t^{1-\alpha} \mbox{ with }0<\alpha<1 \label{eq:prod}
\end{equation}
$a_t$ represents a stochastic shock to technology or Solow residual. We assume that the shocks to technology are distributed with zero mean, but display both persistence across time and correlation in the current period.
Let us consider the joint process $(a_t,b_t)$ defined as
\begin{equation}
\left(
\begin{array}{c}
a_t\\b_t
\end{array}
\right)=
\left(
\begin{array}{cc}
\rho&\tau\\
\tau&\rho\\ 
\end{array}
\right)\left(
\begin{array}{c}
a_{t-1}\\b_{t-1}
\end{array}
\right)
+\left(
\begin{array}{c}
\varepsilon_t\\ \nu_t
\end{array}
\right)  \label{eq:process}
\end{equation}
where $|\rho+\tau|<1$ and $|\rho-\tau|<1 $ for sake of stationarity and 
\begin{eqnarray*}
E(\varepsilon_t)&=& 0,\\
E(\nu_t)&=& 0,\\
E(\varepsilon_t\varepsilon_s)&=&\left\{
\begin{array}{lcl}
\sigma^2_\varepsilon & \mbox{ if } & t=s \\
0 & \mbox{ if } & t\neq s \\
\end{array}\right. \mbox{, }\\ 
E(\nu_t\nu_s)&=&\left\{
\begin{array}{lcl}
\sigma^2_\nu & \mbox{ if } & t=s \\
0 & \mbox{ if } & t\neq s \\
\end{array}\right. \mbox{, } \\
E(\varepsilon_t\nu_s)&=&\left\{
\begin{array}{lcl}
\varphi\sigma_\varepsilon\sigma_\nu & \mbox{ if } & t=s \\
0 & \mbox{ if } & t\neq s \\
\end{array}\right. \mbox{. } 
\end{eqnarray*}

 
\section{Dynamic Equilibrium}
The dynamic equilibrium of this economy follows from the first order conditions for optimality: 
\begin{eqnarray*}
&&c_t \theta h_t^{1+\psi}=(1-\alpha) y_t \\
&&\beta E_t\left[\left(\frac{\exp(b_t) c_t}{\exp(b_{t+1})c_{t+1}}\right)\left(\exp(b_{t+1})\alpha \frac{y_{t+1}}{k_{t+1}}+1-\delta\right)\right]=1\\
&&y_t=\exp(a_t) k_t^\alpha h_t^{1-\alpha} \\
&&k_{t+1}=\exp(b_t)(y_t-c_t)+(1-\delta)k_t \\
&&a_t=\rho a_{t-1}+\tau b_{t-1}+\varepsilon_t \\
&&b_t=\tau a_{t-1}+\rho b_{t-1}+\nu_t 
\end{eqnarray*}
\section{The {\sc dynare} code}
The {\sc dynare} code is straightforward to write, as the equilibrium is written in the natural way. The whole code is reported at the end of the section. Before that we proceed step by step.
\paragraph{Preamble}
The preamble consists of the some declarations to setup the endogenous and exogenous variables, the parameters and assign values to these parameters. 
\begin{enumerate}
\item {\tt var y, c, k, h, a, b;} specifies the endogenous variables in the model since we have output ({\tt y}), consumption ({\tt c}), capital ({\tt k}), hours ({\tt h}) and the two shocks ({\tt a, b}).
\item {\tt varexo e, u;} specifies the exogenous variables in the model --- namely the innovations of the shocks, since we have the innovation of the non--incorporated shock ({\tt e}), and the innovation of the incorporated shock ({\tt u}).
\item {\tt parameters list;} specifies the list of parameters of the model. In the case we are studying:

{\tt parameters beta, alpha, delta, theta, psi, rho, tau;}

\begin{tabular}{ll}
\hline
\hline
{\tt beta}  & discount factor                               \\
{\tt alpha} & capital elasticity in the production function \\
{\tt delta} & depreciation rate                             \\
{\tt theta} & disutility of labor parameter                 \\
{\tt psi}   & labor supply elasticity                       \\
{\tt rho}   & persistence                                   \\
{\tt tau}   & cross--persistence                            \\
\hline
\hline
\end{tabular}
\item Assignment of parameter values. This is done the standard way in MATLAB. For example, we write

\begin{verbatim}
alpha = 0.36;
rho   = 0.95;
tau   = 0.025;
beta  = 0.99;
delta = 0.025;
psi   = 0;
theta = 2.95;
\end{verbatim}


\item Note that $\varphi$, the conditional correlation of the shocks, is not, strickly speaking, a parameter of the recursive equations and doesn't need to be listed in the {\tt parameters} instruction. It may however be convenient to express it as a parameter in the expression of the variance--covariance matrix of the shocks (see below) and one may simply write:

{\tt phi = 0.1; }

\end{enumerate}

\paragraph{Declaration of the model:}
This step is done in a straightforward way. It starts with the instruction {\tt model;} and ends with {\tt end;}, in between all equilibrium conditions are written exactly the way we write it ``by hand''. However, there is a simple rule that should be kept in mind when the model is written. Let us consider a variable $x$: 
\begin{itemize}
\item If $x$ is decided in period $t$ then we simply write ${\tt x}$. 
\item When the variable is decided in $t-1$, such as the capital stock in our simple model, we write ${\tt x(-1)}$. \item Finally, when a variable is decided in the next period, $t+1$, such as consumption in the Euler equation, we write ${\tt x(+1)}$.
\end{itemize}
Hence the required code to declare our model in {\sc Dynare} will be:
\begin{verbatim}
model;
c*theta*h^(1+psi)=(1-alpha)*y;
k = beta*(((exp(b)*c)/(exp(b(+1))*c(+1)))*
    (exp(b(+1))*alpha*y(+1)+(1-delta)*k));
y = exp(a)*(k(-1)^alpha)*(h^(1-alpha));
k = exp(b)*(y-c)+(1-delta)*k(-1);
a = rho*a(-1)+tau*b(-1) + e;
b = tau*a(-1)+rho*b(-1) + u;
end;
\end{verbatim}
Assume now that we want to take a Taylor series expansion in logs rather than in level, we just rewrite the model as
\begin{verbatim}
model;
exp(c)*theta*exp(h)^(1+psi)=(1-alpha)*exp(y);
exp(k) = beta*(((exp(b)*exp(c))/(exp(b(+1))*exp(c(+1))))
         *(exp(b(+1))*alpha*exp(y(+1))+(1-delta)*exp(k)));
exp(y) = exp(a)*(exp(k(-1))^alpha)*(exp(h)^(1-alpha));
exp(k) = exp(b)*(exp(y)-exp(c))+(1-delta)*exp(k(-1));
a = rho*a(-1)+tau*b(-1) + e;
b = tau*a(-1)+rho*b(-1) + u;
end;
\end{verbatim}
so that the level of consumption is actually given by ${\tt exp(c)}$.
\paragraph{Solving the model}
\begin{enumerate} 
\item Now we need to provide numerical initial conditions for the computation of the deterministic steady state. This is done with the sequence between {\tt initval;} and {\tt end;}. Each variable, endogenous or exogenous, should be initialized. In our example, we give the exact values of the deterministic equilibrium in absence of shocks. This takes the form
\begin{verbatim}
initval;
y = 1.08068253095672;
c = 0.80359242014163;
h = 0.29175631001732;
k = 11.08360443260358;
a = 0;
b = 0;
e = 0;
u = 0;
end;
\end{verbatim}

Alternatively, we could provide only approximated values. {\sc Dynare} would then automatically compute the exact values.

\item We then specify the innovations and their matrix of variance--covariance. This is done using a {\tt shocks;} and {\tt end;} block and declare only the nonzero elements of the covariance matrix:
\begin{verbatim}
shocks;
var e; stderr 0.009;
var u; stderr 0.009;
var e, u = phi*0.009*0.009;
end;
\end{verbatim}
Note that it is possible to shut down a shock by assigning it a zero variance.

\item The model is then solved and simulated using the {\tt stoch\_simul;} command. By default, the coefficients of the approximated decision rules are reported as well as the moments of the variables and impulse response functions for each exogenous shocks are ploted. In addition, the following options are available (out of many others):
\begin{itemize}
\item {\tt ar} = Integer:
Order of autocorrelation coefficients to compute and to print (default = 5)
\item {\tt periods} = Integer:
If different from zero, the model will be simulated and empirical moments will be computed instead of theoretical moments. The value of the option specifies the number of periods to use in the simulations (default = 0)
\item {\tt nocorr}:
Doesn't print the correlation matrix (default = PRINT)
\item {\tt drop} = Integer:
Number of points dropped at the beginning of simulation before computing the summary statistics (default = 100)
\item {\tt irf} = Integer:
Number of periods on which to compute the IRFs (default = 40)
\item {\tt nofunctions}:
Doesn't print the coefficients of the approximated solution
\item {\tt nomoments}:
Doesn't print moments of the endogenous variables 
\item {\tt order} = [1,2,3]:
Order of Taylor approximation (default = 2)
\item {\tt replic} = Integer:
Number of simulated series used to compute the IRFs (default = 1, if order = 1, and 50 otherwise)
\end{itemize} 

In our first example, we use simply:
\begin{verbatim}
stoch_simul;
\end{verbatim}
{\sc Dynare} will compute theoretical moments of variables.

In our second example, we use:
\begin{verbatim}
stoch_simul(periods=2000, drop=200);
\end{verbatim}
{\sc Dynare} will compute simulated moments of variables. The simulated trajectories are returned in MATLAB vectors named as the variables (be careful not to use MATLAB reserved names such as INV for your variables \ldots).
\end{enumerate}

\appendix
\section{{\sc Dynare} code for the model in level}
Here is the model file for the model in level. It can be found in file {\tt example1.mod}.

\begin{verbatim}
var y, c, k, a, h, b;
varexo e, u;

parameters beta, rho, alpha, delta, theta, psi, tau;

alpha = 0.36;
rho   = 0.95;
tau   = 0.025;
beta  = 0.99;
delta = 0.025;
psi   = 0;
theta = 2.95;

phi   = 0.1;

model;
c*theta*h^(1+psi)=(1-alpha)*y;
k = beta*(((exp(b)*c)/(exp(b(+1))*c(+1)))
    *(exp(b(+1))*alpha*y(+1)+(1-delta)*k));
y = exp(a)*(k(-1)^alpha)*(h^(1-alpha));
k = exp(b)*(y-c)+(1-delta)*k(-1);
a = rho*a(-1)+tau*b(-1) + e;
b = tau*a(-1)+rho*b(-1) + u;
end;

initval;
y = 1.08068253095672;
c = 0.80359242014163;
h = 0.29175631001732;
k = 11.08360443260358;
a = 0;
b = 0;
e = 0;
u = 0;
end;

shocks;
var e; stderr 0.009;
var u; stderr 0.009;
var e, u = phi*0.009*0.009;
end;

stoch_simul;
\end{verbatim}

\section{{\sc Dynare} code for the model in logs}

Here is the model file for the model in logs. In this case, {\tt initval} only contains guessed values and {\tt steady} is used to compute and display the exact value of the deterministic equilibrium. The shocks are supposed to be uncorrelated. The model file can be found in {\tt example2.mod}.

\begin{verbatim}
var y, c, k, a, h, b;
varexo e, u;

parameters beta, rho, alpha, delta, theta, psi, tau;

alpha = 0.36;
rho   = 0.95;
tau   = 0.025;
beta  = 0.99;
delta = 0.025;
psi   = 0;
theta = 2.95;

model;
exp(c)*theta*exp(h)^(1+psi)=(1-alpha)*exp(y);
exp(k) = beta*(((exp(b)*exp(c))/(exp(b(+1))*exp(c(+1))))
         *(exp(b(+1))*alpha*exp(y(+1))+(1-delta)*exp(k)));
exp(y) = exp(a)*(exp(k(-1))^alpha)*(exp(h)^(1-alpha));
exp(k) = exp(b)*(exp(y)-exp(c))+(1-delta)*exp(k(-1));
a = rho*a(-1)+tau*b(-1) + e;
b = tau*a(-1)+rho*b(-1) + u;
end;

initval;
y = 0.1;
c = -0.2;
h = -1.2;
k =  2.4;
a = 0;
b = 0;
e = 0;
u = 0;
end;

steady;

shocks;
var e = 0.009^2;
var u = 0.009^2;
end;

stoch_simul(periods=2000, drop=200);
\end{verbatim}
\bibliographystyle{Usmad}
\bibliography{/papers/biblio/michel}

\end{document}