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

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

42 43 44 45 | ```
#for i=1:N
var_state z#i;
z#i = linspace(zMin,zMax,zPts); % placeholder
#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.