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

\[\mathbb{E}_0 \sum_{t=0}^{\infty}\beta^t u(c_{j,t}).\]

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

\[Y_{j,t}=\exp(z_{j,t})K_{j,t}^{\alpha},\]

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:

\[K_{j,t+1} =(1-\delta)K_{j,t} + i_{j,t}.\]

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:

\[\begin{split}\max_{c_{j,t},i_{j,t},K_{j,t+1}} \mathbb{E}_0 \sum_j \omega_j \sum_{t=0}^{\infty} \beta^t u(c_{j,t}), \\ s.t. \quad \sum_j (c_{j,t}+i_{j,t})=\sum_j \exp(z_{j,t})K_{j,t}^{\alpha}, \\ K_{j,t+1} =(1-\delta)K_{j,t} + i_{j,t}, \\ \Rightarrow s.t. \quad \sum_j (c_{j,t}+K_{j,t+1})=\sum_j \exp(z_{j,t})K_{j,t}^{\alpha}+(1-\delta)K_{j,t},\end{split}\]

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

\[\begin{split}\omega_j u'(c_{j,t})=\lambda_t, \\ -\lambda_t +\beta E_t \lambda_{t+1}[MPK_{j,t+1} + (1-\delta)]=0,\end{split}\]

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:

\[\begin{split}-\lambda_t + \beta \mathbb{E}_t \lambda_{t+1}[MPK_{j,t+1} + (1-\delta)]=0, \forall j \\ \sum_j (c_{j,t}+K_{j,t+1})=\sum_j \exp(z_{j,t})K_{j,t}^{\alpha}+(1-\delta)K_{j,t},\end{split}\]

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

\[c_{j,t}=u'^{-1}(\frac{\lambda_t}{\omega_j}).\]

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).

% 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(:,GDSGE_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

#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

\[\mathcal{E}_{j,t}=\Big|-1+\frac{ \beta \mathbb{E}_t u'(c_{j,t+1})[MPK_{j,t+1} + (1-\delta)]}{u'(c_{j,t})}\Big|.\]

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.