# Getting Started - A Simple RBC Model¶

## The Model¶

The model is a canonical complete-markets Real Business Cycle (RBC) model with a single productivity shock. This simple model can be approximated with linear solutions well, so using GDSGE to obatin global non-linear solutions is not necessary. Nevertheless it serves as a good starting point to illustrate the basic usage of the toolbox, as this model should be familiar to any macroeconomic audience.

Time is discrete \(t=0...\infty\). Representative households with mass one maximize utility

where parameter \(\beta>0\) is the discount factor, \(\sigma>0\) is the coefficient of relative risk aversion. \(c_t\) is the consumption. \(\mathbb{E}\) is an expectation operator integrating the productivity shock to be introduced below.

Competitive representative firms produce a single output according to

where parameter \(\alpha \in (0,1)\) is the capital share. \(K_t\) is capital, \(L_t\) is labor, and \(z_t\in\{z_L,z_H\}\) is the productivity shock which follows a two-point Markov process with the following transition matrix

Capital depreciates at rate \(\delta \in (0,1)\). Investment technology can convert one unit of output to new capital:

Assets markets are complete. Households make investment decision and supply capital and labor to firms. Market clearing conditions for labor and goods are, respectively,

The complete-markets equilibrium can be characterized by the Euler equation and the goods market clearing condition:

Given \(K_0\), a sequential competitive equilibrium is stochastic processes \((c_t,K_{t+1})\) that satisfy system (1). To input the system into GDSGE, we need to write the system in a recursive form.

Definition: a recursive competitive equilibrium is functions \(c(z,K), K'(z,K)\) s.t.

The recursive system can be solved using a time iteration procedure:

taking function \(c_{t+1}(z,K)\) as known in the period-\(t\) time step.

## The gmod File¶

The recursive system can now be input to GDSGE via a mod file `rbc.gmod`

below.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 | ```
% Parameters
parameters beta sigma alpha delta;
beta = 0.99; % discount factor
sigma = 2.0; % CRRA coefficient
alpha = 0.36; % capital share
delta = 0.025; % depreciation rate
% Exogenous States
var_shock z;
shock_num = 2;
z_low = 0.99; z_high = 1.01;
Pr_ll = 0.9; Pr_hh = 0.9;
z = [z_low,z_high];
shock_trans = [
Pr_ll, 1-Pr_ll
1-Pr_hh, Pr_hh
];
% Endogenous States
var_state K;
Kss = (alpha/(1/beta - 1 + delta))^(1/(1-alpha));
KPts = 101;
KMin = Kss*0.9;
KMax = Kss*1.1;
K = linspace(KMin,KMax,KPts);
% Interp
var_interp c_interp;
initial c_interp z.*K.^alpha+(1-delta)*K;
% Time iterations update
c_interp = c;
% Endogenous variables as unknowns of equations
var_policy c K_next;
inbound c 0 z.*K.^alpha+(1-delta)*K;
inbound K_next 0 z.*K.^alpha+(1-delta)*K;
% Other endogenous variables
var_aux w;
model;
% Budget constraints
u_prime = c^(-sigma);
kret_next' = z'*alpha*K_next^(alpha-1) + 1-delta;
% Evaluate the interpolation object to get future consumption
c_future' = c_interp'(K_next);
u_prime_future' = c_future'^(-sigma);
% Calculate residual of the equation
euler_residual = 1 - beta*GDSGE_EXPECT{u_prime_future'*kret_next'}/u_prime;
market_clear = z*K^alpha + (1-delta)*K - c - K_next;
% Calcualte other endogenous variables
w = z*(1-alpha)*K^alpha;
equations;
euler_residual;
market_clear;
end;
end;
simulate;
num_periods = 10000;
num_samples = 6;
initial K Kss;
initial shock 1;
var_simu c K w;
K' = K_next;
end;
``` |

The gmod file can be uploaded to the online compiler at http://www.gdsge.com/ (Remember also to download the necessary runtime files following the instruction on the compiler website). The compiler returns three files that can be used to solve and simulate the model: iter_rbc.m, simulate_rbc.m, mex_rbc.mexw64.

First, call iter_rbc.m in matlab to run the policy iterations and store the returned result in *IterRslt*, which produces

```
>> IterRslt = iter_rbc;
Iter:10, Metric:0.385607, maxF:9.99011e-09
Elapsed time is 0.174117 seconds.
...
Iter:323, Metric:9.89183e-07, maxF:7.96886e-09
Elapsed time is 0.027923 seconds.
```

where *Metric* is the maximum absolute distance of *var_interp* between the current and the last iteration,
and *maxF* is the maximum absolute equation residual across all collocation points in the current iteration.
As shown, the converged criterion by default is *Metric* smaller than \(1e-6\), and can be overwritten by the option *TolEq* that can be put into the gmod file or
supplied in an option structure when calling *iter_rbc* (see details below).

The returned IterRslt contains the converged policy and state transition functions. For example, we can plot the policy function for consumption \(c\):

```
>> plot(IterRslt.var_state.K, IterRslt.var_policy.c);
xlabel('K'); title('Policy Functions for c');
```

We can now simulate the model by inputting IterRslt into simulate_rbc:

```
>> SimuRslt = simulate_rbc(IterRslt);
Periods: 1000
shock K c w
2 36.27 2.707 2.448
Elapsed time is 0.642607 seconds.
...
Periods: 10000
shock K c w
1 36.71 2.685 2.224
Elapsed time is 0.632642 seconds.
```

The returned SimuRslt contains a panel of simulated paths with num_samples and num_periods specified in the mod file. For example, we can plot the histogram of state variable \(K\):

```
>> histogram(SimuRslt.K); title('Histogram for K');
```

Or we can plot the first two sample paths of wages for the first 100 periods:

```
>> plot(SimuRslt.w(1:2,1:100)'); title('Sample Paths of Wages');
```

The iter and simulate files can be reused by passing parameters to be overwritten in a struct. For example,

```
>> options.sigma = 1.5 % overwrite sigma
options.z = [0.95,1.05] % making the shock larger
IterRslt = iter_rbc(options);
SimuRslt = simulate_rbc(IterRslt,options);
```

Part of the toolbox options can also be overwritten by including them in the struct. See Options.

## Dissecting the gmod File¶

We now look into each line of rbc.gmod and describe the usage.

1 2 3 4 5 6 | ```
% Parameters
parameters beta sigma alpha delta;
beta = 0.99; % discount factor
sigma = 2.0; % CRRA coefficient
alpha = 0.36; % capital share
delta = 0.025; % depreciation rate
``` |

A line starting with keyword parameters declares parameters.
Parameters are variables used in the *model* block that are invariant across states or over time.
Matlab-style inline comments can be used following a “%”.

8 9 10 11 12 13 14 15 16 17 | ```
% Exogenous States
var_shock z;
shock_num = 2;
z_low = 0.99; z_high = 1.01;
Pr_ll = 0.9; Pr_hh = 0.9;
z = [z_low,z_high];
shock_trans = [
Pr_ll, 1-Pr_ll
1-Pr_hh, Pr_hh
];
``` |

Exogenous states such as the productivity shock in the RBC model should be declared following keyword *var_shock*.
The number of states should be specified in variable *shock_num*.
The Markov transition matrix should be specified in variable *shock_trans*.
Notice that any Matlab expressions such as `z_low = 0.99`

can be evaluated outside the *model* block.

For multi-dimension exogenous states, *shock_num* needs to be specified as the number of discretze points in the Cartesian set across
all dimensions. *shock_trans* should be specified as the full transition matrix accordingly.
Matlab functions *ndgrid* and *kron* can be used to construct the Cartesian set and the full transition matrix.

If an exogenous state follows a continuous process such as AR(1), specify the state to be an endogenous state instead in *var_state*,
and discretize the innovation.

19 20 21 22 23 24 25 | ```
% Endogenous States
var_state K;
Kss = (alpha/(1/beta - 1 + delta))^(1/(1-alpha));
KPts = 101;
KMin = Kss*0.9;
KMax = Kss*1.1;
K = linspace(KMin,KMax,KPts);
``` |

Endogenous states such as capital in the RBC model should be declared following keyword *var_state*.
Then the discretized grid for each state should be specified (e.g., `K = linspace(KMin,KMax,KPts);`

in the example).
Notice that any Matlab expressions (such as function *linspace* here) can be evaluated.

The discretized grid will be used in fixed-grid function approximation procedures such as splines and linear interpolations. For adaptive grid methods, only the range of the grid will be used.

27 28 29 30 31 | ```
% Interp
var_interp c_interp;
initial c_interp z.*K.^alpha+(1-delta)*K;
% Time iterations update
c_interp = c;
``` |

The time iteration procedure requires taking the state transition functions (here, implicitly defined in \(c_{t+1}(z,K)\)) as given
when evaluating the system of equations. These state transition functions are declared following keyword *var_interp*.
They are named *var_interp* since they are usually evaluated with some function approximation procedure such as interpolation in evaluating the system of equations.

Each of the transition functions needs to be initialized following keyword *initial*.
Here \(c_{t+1}(z,K)\) is initialized as \(zK^{\alpha} + (1-\delta)K\) which means that consumers just
consumer all resources that include output and non-depreciated capital. This initialization basically
requests the toolbox to find the equilibrium that corresponds to the limit of a finite-horizon economy taking the number of horizons to infinity,
of which \(c_{t+1}\) is initialized to be the last-period solution.
Notice that the Matlab “dot” operator (.*) works in the line following *initial*, and relevant tensors are formed automatically.

In general, it is crucial to form a good initial guess on the transition functions to make the policy iteration method work.
Starting with a last-period problem is shown to
deliver both good theoretical properties and robust numerical computations (Cao, 2018).
See Mendoza (2010) for an example on how to define a more complex last-period problem that may require solving a different system of equations, through a *model_init* block.

The update of the transition function after a time iteration step needs to be specified such as `c_interp = c;`

.
The update can use any variables returned as part of the solution to the system in the current time iteration step.
For example, here `c_interp = c;`

means that the *c_interp* variable is updated with *c* after a time iteration,
where *c* is one of the policy variables solved by the system of equations, specified in the gmod file as below.

33 34 35 36 | ```
% Endogenous variables as unknowns of equations
var_policy c K_next;
inbound c 0 z.*K.^alpha+(1-delta)*K;
inbound K_next 0 z.*K.^alpha+(1-delta)*K;
``` |

All policy variables that enter the equations as unknowns are declared following keyword `var_policy`

.
For each of these variables, the bounds for where the equation solver searches the solution should be specified
following keyword *inbound*. For variables of which the bounds cannot be determined ex-ante,
specify tight bounds for them and use an *adaptive* option so the solver adjusts bounds automatically when it cannot find solutions
and adapt bounds across time iterations. See `var_policy`

in Variable declaration.

38 39 | ```
% Other endogenous variables
var_aux w;
``` |

Some policy variables can be evaluated as simple functions of states and other policy variables, so they do not have to be
directly included in the system of equations. For example wage is determined by the marginal product of labor
and can be evaluated according to \(w=(1-\alpha)zK^{\alpha}\). These policy variables can be declared
following keyword `var_aux`

.

41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 | ```
model;
% Budget constraints
u_prime = c^(-sigma);
kret_next' = z'*alpha*K_next^(alpha-1) + 1-delta;
% Evaluate the interpolation object to get future consumption
c_future' = c_interp'(K_next);
u_prime_future' = c_future'^(-sigma);
% Calculate residual of the equation
euler_residual = 1 - beta*GDSGE_EXPECT{u_prime_future'*kret_next'}/u_prime;
market_clear = z*K^alpha + (1-delta)*K - c - K_next;
% Calcualte other endogenous variables
w = z*(1-alpha)*K^alpha;
equations;
euler_residual;
market_clear;
end;
end;
``` |

The *model;—end;* block defines the system of equations for each collocation point of endogenous states and exogenous states.
For the current example, it simply defines the system of equations for each collocation point \((z,K)\).

The equations should be eventually specified in the *equation;—end;* block, in which each line corresponds to one equation in the system.
Any calculations in order to evaluate these equations are included in the *model* block preceding the *equations* block.
Notice the whole *model* block is parsed into C++, so all evaluations should be scalar-based: Matlab functions such as .* operator cannot be used.
Nevertheless, the *model* block supports simple control-flow blocks in Matlab-style syntax, such as if-elseif-else-end block and
for-end block.

A variable followed by a prime (‘) indicates that the variable is a vector of length *shock_num*.
For example,

44 | ```
kret_next' = z'*alpha*K_next^(alpha-1) + 1-delta;
``` |

means that

```
kret_next(1) = z(1)*alpha*K_next^(alpha-1) + 1-delta
kret_next(2) = z(2)*alpha*K_next^(alpha-1) + 1-delta
```

A var_shock (here `z`

), when used followed by a prime (‘), corresponds to the vector of the variable across all possible exogenous states, and thus can be
used to construct other vector variables. Notice when a var_shock is used not followed by a prime, it corresponds to the exogenous state at the current collocation point, e.g., `z`

in

52 | ```
market_clear = z*K^alpha + (1-delta)*K - c - K_next;
``` |

Any *var_interp* declared before can be used as a function in the model block. For example,

47 | ```
c_future' = c_interp'(K_next);
``` |

means calculating \(c_{t+1}\) based on the state transition function implied by *c_interp*, at
the next-period endogenous state *K_next* (which is a current policy variable).
Where are the next-period exogenous states? These are taken care of by the prime (‘) operator
following *c_interp*, which means that it is evaluating at each realization of the future exogenous states,
and returning a vector of length shock_num. Accordingly, the left hand side variable should be declared as a vector,
i.e., a variable followed by a prime (‘).

The *model* block can use several built-in functions for reduction operations.
For example, the *GDSGE_EXPECT{}* used in

51 | ```
euler_residual = 1 - beta*GDSGE_EXPECT{u_prime_future'*kret_next'}/u_prime;
``` |

is to calculate the conditional expectation of the object defined inside the curly braces,
conditional on the current exogenous state, using the default Markov transition matrix *shock_trans*.
Obviously, this function is meaningful only if it takes as argument the realizations of variables across future exogenous states,
which are defined as vector variables followed by prime (‘).

This operator can also take a different transition matrix than *shock_trans*, which is used as
`GDSGE_EXPECT{*|alternative_shock_trans}`

. This can be used to solve models with heterogeneous beliefs.
See example Cao (2018).

Two other reduction operations *GDSGE_MAX{}* and *GDSGE_MIN{}* are defined, which are to take
the maximum and the minimum of objects inside the curly braces, respectively. See Built-in functions.

51 52 53 54 55 56 57 58 | ```
simulate;
num_periods = 10000;
num_samples = 6;
initial K Kss;
initial shock 1;
var_simu c K w;
K' = K_next;
end;
``` |

The simulate block defines how the file for Monte Carlo simulations should be generated.
It should define the initial endogenous state following keyword *initial*.
`initial K Kss;`

in the example sets the endogenous state \(K\) to *Kss* defined before.

The simulate block should define the initial exogenous state index following keyword *initial shock*.
It should define the variables to be recorded following *var_simu*.
It should define the transition for each state variable. In the example, `K'=K_next;`

defines that the next period endogenous state \(K\) should be assigned to *K_next* which is one of
the *var_policy* solved as part of the system. Notice the prime operator (‘) following *K*
indicates that the line is to specify the transition of an endogenous state variable in the simulation,
which has a different meaning than the prime operator (‘) used in the *model* block.
Nevertheless, the prime operators (‘) in both usages are associated with the transition to the future state, and thus motivate such designs.

The simulate block can also overwrite num_periods (default 1000) and num_samples (default 1).

## What’s Next?¶

The current example demonstrates the basic usage of the toolbox. You can proceed to an extension with investment irreversibility that requires a global method, or to a real example Heaton and Lucas (1996) which is the leading example in the toolbox paper.

Or you can proceed to Toolbox API.