# A Multi-country Business Cycle Model

## The Model

In this example we demonstrate the toolbox’s ability to handle models with high dimensional states, thanks to the sparse-grid interpolation method. The model is a complete-markets multi-country business cycle model. We show that the model with 15 countries and 30 continuous states can be solved in 10 minutes with decent accuracy on a regular desktop.

There are \(N\) countries, indexed by \(j\). Each country has one unit of consumers that supply labor exogenously. The consumers’ preference reads

Production in a country \(j\) takes capital and labor according to a Cobb-Douglas production function:

where \(K_{j,t}\) is capital and \(z_{j,t}\) is the log of country-j’s productivity that follows an AR(1) process, which is independently distributed across countries. Capital and labor are immobile. Capital depreciates at rate \(\delta\) and investment follows a standard technology:

The market is complete in that consumers can trade a full menu of Arrow securities.

## Equilibrium Conditions

The complete-markets equilibrium can be characterized by the social planner’s solution:

where \(\omega_j\) is the Pareto weight assigned to consumers in country \(j\) that depends on the initial distribution of Arrow securities. The first order conditions read

where \(MPK_{j,t}\equiv \alpha \exp(z_{j,t})K_{j,t}^{\alpha-1}\), and \(\lambda_t\) is the Lagrangian multiplier associated with the world resource constraint. Therefore, the equilibrium conditions can be stated as:

for equilibrium objects \((\lambda_t,K_{j,t+1})\), where

## Solutions

We solve the recursive equilibrium defined over the state vector \((z_1,K_1,z_2,K_2,...,z_N,K_N)\) of length \(2\times N\). The model does not admit an analytical solution despite being simple. The numerical solution is challenging as the dimension of state space is high, and interpolation methods based on dense grids face the curse of dimensionality. In particular, the number of collocation points increases exponentially in the dimension if dense-grid methods are used. Instead, we use the adaptive-grid method that is shipped with the toolbox, with which the number of collocation points only increases in the dimension in a polynomial speed.

The gmod code is listed below where one can find the detailed parameterization (`rbc`

).

```
1% Number of countries
2#define N 15
3
4% Use asg
5USE_SPLINE = 0; USE_ASG = 1;
6AsgMinLevel = 2;
7AsgMaxLevel = 2;
8AsgThreshold = 1e-2;
9AsgOutputMaxLevel = 4;
10AsgOutputThreshold = 1e-2;
11PrintFreq = 1;
12SaveFreq = inf;
13TolEq = 1e-5;
14SIMU_RESOLVE = 0; SIMU_INTERP = 1;
15
16% Parameters
17parameters beta sigma alpha delta numCountries;
18beta = 0.99; % discount factor
19sigma = 2.0; % CRRA coefficient
20alpha = 0.36; % capital share
21delta = 0.025; % depreciation rate
22numCountries = N;
23
24% Innovations
25parameters rho_z sigma_z prob_z;
26rho_z = 0.9; % First-order autocorrelation coefficient
27sigma_z = 0.01; % std of innovation of productivity
28prob_z = 1/(2*N); % Probability for each perturbation; used for integration
29
30% Pareto weights
31parameters omega_vec;
32omega_vec = ones(1,N);
33
34% Endogenous States
35Kss = (alpha/(1/beta - 1 + delta))^(1/(1-alpha));
36KPts = 2; % placeholder
37KMin = Kss*0.9;
38KMax = Kss*1.1;
39zMin = -0.05;
40zMax = 0.05;
41zPts = 2; % placeholder
42#for i=1:N
43var_state z#i;
44z#i = linspace(zMin,zMax,zPts); % placeholder
45#end
46#for i=1:N
47var_state K#i;
48K#i = linspace(KMin,KMax,KPts); % placeholder
49#end
50
51var_policy_init lambda;
52inbound_init lambda 0 1;
53model_init;
54 % Construct vector from state
55 vector z_vec[N];
56 vector K_vec[N];
57 #for i=1:N
58 z_vec(#i) = z#i;
59 K_vec(#i) = K#i;
60 #end
61
62 % Back out consumption from the marginal utility
63 vector c_vec[N];
64 for i=1:N
65 c_vec(i) = (lambda / omega_vec(i))^(-1/sigma);
66 end
67
68 % Calcualte resource residual
69 resource_resid = 0;
70 for i=1:N
71 resource_resid = resource_resid + exp(z_vec(i))*K_vec(i)^alpha + (1-delta)*K_vec(i) - c_vec(i) - K_vec(i);
72 end
73 equations;
74 resource_resid;
75 end;
76end;
77
78% Interp
79var_interp lambda_interp;
80initial lambda_interp lambda;
81lambda_interp = lambda;
82
83% Endogenous variables as unknowns of equations
84var_policy K_next[N]
85inbound K_next 25 45;
86var_policy lambda;
87inbound lambda 0 1;
88
89model;
90 % Construct vector from state
91 vector z_vec[N];
92 vector K_vec[N];
93 #for i=1:N
94 z_vec(#i) = z#i;
95 K_vec(#i) = K#i;
96 #end
97
98 % Prepare future states, under no perturbation
99 vector state_next[2*N];
100 for i=1:N
101 % First N states is future productivity
102 state_next(i) = rho_z*z_vec(i);
103 % Next N states if future capital
104 state_next(N+i) = K_next(i);
105 end
106
107 % Calculate marginal product of capital with no perturbation
108 vector mpk_next[N];
109 for i=1:N
110 mpk_next(i) = alpha*exp(state_next(i)) * K_next(i)^(alpha-1);
111 end
112
113 % Do a sparse quadrature by perturbing productivity once at a time
114 vector euler_residual[N];
115 for i=1:N
116 euler_residual(i) = 0;
117 end
118 for i=1:N
119 % Modifty future productivity of the perturbed country
120 % Positive innovation
121 state_next(i) = rho_z*z_vec(i) + sigma_z*sqrt(numCountries);
122 mpk_next(i) = alpha*exp(state_next(i)) * K_next(i)^(alpha-1);
123 lambda_future = GDSGE_INTERP_ARRAY_adouble(shock,state_next_GDSGE_GRID);
124 #for j=1:N
125 euler_residual(#j) = euler_residual(#j) + prob_z*(-lambda + beta*lambda_future*(mpk_next(#j)+(1-delta)));
126 #end
127
128 % Negative innovation
129 state_next(i) = rho_z*z_vec(i) - sigma_z*sqrt(numCountries);
130 mpk_next(i) = alpha*exp(state_next(i)) * K_next(i)^(alpha-1);
131 lambda_future = GDSGE_INTERP_ARRAY_adouble(shock,state_next_GDSGE_GRID);
132 #for j=1:N
133 euler_residual(#j) = euler_residual(#j) + prob_z*(-lambda + beta*lambda_future*(mpk_next(#j)+(1-delta)));
134 #end
135
136 % Revert the perturbation back
137 state_next(i) = rho_z*z_vec(i);
138 mpk_next(i) = alpha*exp(state_next(i)) * K_next(i)^(alpha-1);
139 end
140
141 % Back out consumption from the marginal utility
142 vector c_vec[N];
143 for i=1:N
144 c_vec(i) = (lambda / omega_vec(i))^(-1/sigma);
145 end
146 % Calcualte resource residual
147 resource_resid = 0;
148 for i=1:N
149 resource_resid = resource_resid + exp(z_vec(i))*K_vec(i)^alpha + (1-delta)*K_vec(i) - c_vec(i) - K_next(i);
150 end
151
152 % Record variables
153 #for i=1:N
154 K_next_output#i = K_next(#i);
155 c#i = c_vec(#i);
156 #end
157
158 equations;
159 for i=1:N
160 euler_residual(i);
161 end
162 resource_resid;
163 end;
164end;
165
166simulate;
167 num_samples = 24;
168 num_periods = 1000;
169 #for i=1:N
170 initial shock 1;
171 initial z#i 0;
172 initial K#i Kss;
173 var_simu K_next_output#i c#i;
174 z#i' = z_exo(:,GDSGE_t+1,#i);
175 K#i' = K_next_output#i;
176 #end
177end;
```

The code utilizes the macro feature of the toolbox. For example, the code block

```
42#for i=1:N
43var_state z#i;
44z#i = linspace(zMin,zMax,zPts); % placeholder
45#end
```

expands to

```
var_state z1;
z1 = linspace(zMin,zMax,zPts); % placeholder
var_state z2;
z2 = linspace(zMin,zMax,zPts); % placeholder
...
var_state z15;
z15 = linspace(zMin,zMax,zPts); % placeholder
```

The strategy we use for integrating over future shocks follows the monomial integration proposed by Judd, Maliar and Maliar (2011, QE), which calculates the integration by perturbing future innovation from zero along one dimension at a time. This is implemented by Line 118-139 of the code.

The policy functions can be solved by simply calling *iterRslt = iter_rbc.m* after the gmod is compiled. The policy iteration takes around 10 minutes on a regular desktop (with a 12-core 2.5GHz CPU).

Simulation, however, requires some manual configuration as now we need to feed in shocks generated from the actual AR(1) processes, instead of the realizations of shocks used for integration in solving the policy functions. Please refer to `main_simulate`

for how to generate a Markov process using MATLAB’s built-in functions, and feed in the simulated processes into
the toolbox-generated code.

The file also contains procedures in calculating the Euler equation errors. The unit-free Euler equation error is defined as

The error given one state can be calculated by simulating a large sample one period forward starting from this state, and integrating over these sample paths. The script file implements this procedure and reports the max and mean errors for a sample of states drawn from the ergodic distribution (running the calculation requires 16GB memory; set the number of paths lower to reduce the memory requirement). The max error is 6.4E-4 and the mean error is 5.7E-4 across states in the ergodic set, demonstrating that the solutions are obtained with decent accuracy, despite the very sparse grid used for interpolation.

## What’s next

The current example demonstrates that the toolbox can handle models with high dimensional states, with dimension up to 30 on a regular desktop/laptop. However, we note that providing an accurate global solution to a highly-nonlinear albeit low-dimensional model is often much harder than providing a solution to a high-dimensional but less nonlinear model. Nevertheless, this means that one can solve standard medium New-Keynesian models with the toolbox (e.g., the benchmark model in Smets and Wouters (2007) has 9 exogenous and 11 endogenous state variables), and extend these models with nonlinear features such as an occasionally binding interest rate zero lower bound as in the ZLB example.