Government Spending RBC Model for Italy

This notebook calibrates a real business cycle model with government spending using Italian Penn World Table data for 2010 $-$ 2023. The model is then solved in Dynare and compared with empirical business-cycle moments.

1. Model setup

We start from a standard real business cycle model with government spending. Households choose consumption, labor supply and savings, firms produce output using capital and labor, and the government finances exogenous public spending through a labor income tax.

The household Euler equation is

$$ \frac{C_{t+1}}{C_t} = \beta (1+r_t), $$

where $C_t$ is consumption, $\beta$ is the discount factor, and $r_t$ is the real interest rate.

The labor supply condition is

$$ w_t = \frac{\theta}{1-\tau_t}\frac{C_t}{1-N_t}, $$

where $w_t$ is the real wage, $N_t$ is labor input, $\tau_t$ is the labor income tax rate, and $\theta$ controls the weight of leisure in household preferences.

Firms produce output using a Cobb--Douglas technology:

$$ Y_t = A_t K_t^{\alpha} N_t^{1-\alpha}, $$

where $Y_t$ is output, $A_t$ is total factor productivity, $K_t$ is capital, and $\alpha$ is the capital share.

Factor prices are given by marginal products:

$$ w_t = \frac{(1-\alpha)Y_t}{N_t}, $$

$$ r_t+\delta = \frac{\alpha Y_{t+1}}{K_{t+1}}, $$

where $\delta$ is the depreciation rate.

Government spending is financed by labor income taxation:

$$ \tau_t = \frac{G_t}{w_t N_t}, $$

where $G_t$ denotes government spending.

Capital evolves according to

$$ K_{t+1} = (1-\delta)K_t + I_t, $$

and the aggregate resource constraint is

$$ Y_t = C_t + I_t + G_t. $$

Government spending follows an exogenous AR(1) process:

$$ \log(G_t) = (1-\rho_g)\log(\mu^G Y) + \rho_g \log(G_{t-1}) + \varepsilon_{g,t}. $$

Technology also follows an AR(1) process:

$$ \log(A_t) = \rho \log(A_{t-1}) + \varepsilon_t. $$

Main variables

Symbol Meaning
$Y_t$ output
$C_t$ consumption
$I_t$ investment
$G_t$ government spending
$K_t$ capital stock
$N_t$ labor input
$A_t$ total factor productivity
$w_t$ real wage
$r_t$ real interest rate
$\tau_t$ labor income tax rate

Main parameters

Parameter Meaning
$\beta$ discount factor
$\alpha$ capital share
$\delta$ depreciation rate
$\theta$ preference weight on leisure
$\rho$ technology persistence
$\rho_g$ government spending persistence
$\mu^G$ steady-state government spending-to-output ratio

2. Data construction

The next block reads the Penn World Table data, keeps Italy for 2010 $-$ 2023, and constructs the variables used in the model. Output, capital, consumption, investment and government spending are measured in real levels. Labor input is total annual hours worked, and total factor productivity is taken directly from the PWT TFP series.

MATLAB codeCell 1
Code
%% =========================================================
%  PWT 11.0 - Italy, 2010-2023
%  Extract: A, Y, N, K, I, C, G, alpha
%% =========================================================

clear; clc;

%% 1. File path

file_path = "C:\Users\User\Documents\Dynamic_macro\RBC project\pwt110.xlsx";   

%% 2. Read PWT data sheet

T = readtable(file_path, ...
    "Sheet", "Data", ...
    "VariableNamingRule", "preserve");

%% 3. Filter Italy and years 2010-2023

idx = strcmp(T.countrycode, "ITA") & T.year >= 2010 & T.year <= 2023;

IT = T(idx, :);

%% 4. Extract variables

year = IT.year;

% Output
Y = IT.rgdpo;        % Real GDP, output-side, chained PPP

% Capital
K = IT.rnna;         % Capital stock at constant national prices

% Shares of GDP
C_Y = IT.csh_c;      % Consumption share of GDP
I_Y = IT.csh_i;      % Investment share of GDP
G_Y = IT.csh_g;      % Government spending share of GDP

% Levels from shares
C = C_Y .* Y;
I = I_Y .* Y;
G = G_Y .* Y;

% Labor input
emp = IT.emp;        % Persons engaged, millions
avh = IT.avh;        % Average annual hours worked

N = emp .* avh;      % Total annual hours worked, in million hours

% Technology
A = IT.rtfpna;       % TFP at constant national prices

% Capital share
labor_share = IT.labsh;
alpha = 1 - labor_share;

%% 5. Create final table

italy_rbc = table( ...
    year, A, Y, N, K, I, C, G, alpha, ...
    C_Y, I_Y, G_Y, emp, avh, labor_share ...
);

%% 6. Display result

disp(italy_rbc);
Output
Output
    year       A           Y           N          K             I             C             G          alpha       C_Y        I_Y        G_Y       emp       avh      labor_share
    ____    _______    __________    _____    __________    __________    __________    __________    _______    _______    _______    _______    ______    ______    ___________

    2010     1.0566    2.5179e+06    43103    1.8346e+07    6.5211e+05    1.5414e+06    4.3581e+05    0.45962    0.61219    0.25899    0.17309    24.726    1743.2      0.54038  
    2011     1.0557    2.6235e+06    43100    1.8538e+07    6.8081e+05    1.6095e+06    4.3812e+05     0.4636     0.6135    0.25951      0.167    24.784      1739       0.5364  
    2012     1.0315    2.5626e+06    42042    1.8641e+07    6.2105e+05    1.5517e+06    4.2235e+05    0.45597     0.6055    0.24235    0.16481    24.721    1700.6      0.54403  
    2013     1.0239     2.481e+06    40957     1.869e+07    5.7812e+05    1.4858e+06    4.1194e+05      0.462    0.59888    0.23302    0.16603    24.271    1687.5        0.538  
    2014     1.0206     2.476e+06    40914    1.8728e+07    5.8181e+05    1.4592e+06    4.1085e+05    0.46528    0.58935    0.23498    0.16593    24.273    1685.6      0.53472  
    2015     1.0201     2.486e+06    41244    1.8784e+07    5.9447e+05    1.4674e+06    4.0944e+05    0.46369    0.59026    0.23913     0.1647    24.433    1688.1      0.53631  
    2016      1.017    2.6198e+06    41947    1.8868e+07    6.4554e+05    1.5132e+06    4.2128e+05     0.4725    0.57761    0.24641    0.16081    24.789    1692.1       0.5275  
    2017     1.0199    2.6974e+06    42422    1.8969e+07    6.8138e+05    1.5496e+06     4.314e+05    0.46986    0.57447     0.2526    0.15993    25.092    1690.7      0.53014  
    2018      1.015    2.6985e+06    42851     1.909e+07    7.0037e+05    1.5497e+06    4.2308e+05     0.4651    0.57428    0.25954    0.15678    25.352    1690.2       0.5349  
    2019     1.0114    2.8233e+06    42867    1.9216e+07    6.9189e+05    1.6372e+06    4.3894e+05    0.46283     0.5799    0.24506    0.15547    25.509    1680.5      0.53717  
    2020    0.98095    2.5837e+06    37875    1.9276e+07    6.1976e+05    1.4654e+06    4.5043e+05    0.46145    0.56719    0.23988    0.17434    25.018    1513.9      0.53855  
    2021          1    2.8497e+06    41841    1.9507e+07    8.3218e+05    1.5256e+06    4.6309e+05    0.46907    0.53535    0.29202     0.1625    25.312      1653      0.53093  
    2022     1.0083    2.9777e+06    43637    1.9798e+07    9.3839e+05     1.598e+06    4.6386e+05    0.48123    0.53666    0.31514    0.15578    25.802    1691.2      0.51877  
    2023    0.98692     3.004e+06    44672    2.0148e+07    8.9818e+05    1.6193e+06    4.7316e+05    0.48835    0.53903    0.29899    0.15751    26.256    1701.4      0.51165  


The extracted data give the empirical basis for the calibration. The sample is short, so the later business-cycle statistics should be interpreted carefully, but it is enough to obtain Italy-specific average ratios such as consumption, investment, government spending and capital relative to output.

3. Log series

The next figure plots the log levels of the main macroeconomic variables. This is mainly a diagnostic step: it shows the broad movements in productivity, output, capital, consumption, investment and government spending before the calibration is computed.

MATLAB codeCell 2
Code
%% Plot selected RBC variables for Italy

vars_to_plot = {'A','C','I','K','G','Y'};
titles_ = {'TFP (A)', 'Consumption (C)', 'Investment (I)', ...
           'Capital Stock (K)', 'Government Spending (G)', 'Output (Y)'};

%% Plot log series as well

figure;
tiledlayout(3,2, "Padding","compact", "TileSpacing","compact");

for j = 1:numel(vars_to_plot)
    nexttile;
    plot(italy_rbc.year, log(italy_rbc.(vars_to_plot{j})), '-o', 'LineWidth', 1.5);
    grid on;
    xlabel('Year');
    ylabel(['log(' vars_to_plot{j} ')']);
    title(['Log ' titles_{j}]);
end
Output
Output
Output image from Cell 2

The plotted series show the key features of the Italian sample: output and expenditure components fluctuate visibly around the crisis and recovery years, while capital moves more smoothly. This supports using detrended log variables when comparing empirical and model-implied business-cycle moments.

4. Calibration from data ratios

This block computes average ratios from the Italian data: consumption-to-output, investment-to-output, government spending-to-output and capital-to-output. It also backs out the depreciation rate, the capital share, the steady-state real interest rate and the discount factor.

MATLAB codeCell 3
Code
%% =========================================================
%  Calibration ratios from Italian PWT data, 2010-2023
%% =========================================================

% Basic ratios
C_Y = italy_rbc.C ./ italy_rbc.Y;
I_Y = italy_rbc.I ./ italy_rbc.Y;
G_Y = italy_rbc.G ./ italy_rbc.Y;
K_Y = italy_rbc.K ./ italy_rbc.Y;

% Depreciation implied by steady-state relation: I = delta * K
delta_series = italy_rbc.I ./ italy_rbc.K;

% Capital share
alpha_series = italy_rbc.alpha;

%% Averages

c_y_bar     = mean(C_Y, "omitnan");
i_y_bar     = mean(I_Y, "omitnan");
g_y_bar     = mean(G_Y, "omitnan");
k_y_bar     = mean(K_Y, "omitnan");
delta_bar   = mean(delta_series, "omitnan");
alpha_bar   = mean(alpha_series, "omitnan");

%% Implied interest rate and beta

% From RBC steady state:
% r + delta = alpha * Y / K
r_bar = alpha_bar / k_y_bar - delta_bar;

% Euler equation:
% 1 = beta * (1 + r)
beta_bar = 1 / (1 + r_bar);

%% Calibration table

calibration_table = table( ...
    ["C/Y"; "I/Y"; "G/Y"; "K/Y"; "delta"; "alpha"; "r"; "beta"], ...
    [c_y_bar; i_y_bar; g_y_bar; k_y_bar; delta_bar; alpha_bar; r_bar; beta_bar], ...
    'VariableNames', {'Parameter', 'Value'} ...
);

disp(calibration_table);
Output
Output
    Parameter     Value  
    _________    ________

     "C/Y"        0.57816
     "I/Y"         0.2584
     "G/Y"        0.16319
     "K/Y"         7.1469
     "delta"     0.036346
     "alpha"      0.46718
     "r"         0.029022
     "beta"        0.9718


The calibration implies an average consumption share of about 0.58, investment share of about 0.26 and government spending share of about 0.16. The capital-output ratio is high, around 7.15, while the implied annual real interest rate is about 2.9 percent and the discount factor is about 0.972.

5. Shock processes

Technology and government spending are modeled as exogenous AR(1) processes. Technology is estimated from demeaned log TFP. Government spending is first HP-filtered, and the cyclical component is then fitted with an AR(1) process.

MATLAB codeCell 4
Code
%% =========================================================
%  Estimate AR(1) parameters for A and detrended G
%  Italy, 2010-2023
%% =========================================================

lambda = 100;   % HP filter parameter for annual data

%% 1. Technology process: log(A_t)

logA = log(italy_rbc.A);

% Demean log A, so the process is centered around zero
a = logA - mean(logA, "omitnan");

a_lag = a(1:end-1);
a_now = a(2:end);

valid_A = ~isnan(a_lag) & ~isnan(a_now);

rho_A = a_lag(valid_A) \ a_now(valid_A);

eps_A = a_now(valid_A) - rho_A * a_lag(valid_A);
sigma_A = std(eps_A, "omitnan");

%% 2. Government spending process: detrended log(G_t)

logG = log(italy_rbc.G);

% HP filter: logG = trend + cycle
[logG_trend, g_cycle] = hpfilter(logG, "Smoothing", lambda);

% Detrended government spending
g = g_cycle;

g_lag = g(1:end-1);
g_now = g(2:end);

valid_G = ~isnan(g_lag) & ~isnan(g_now);

rho_G = g_lag(valid_G) \ g_now(valid_G);

eps_G = g_now(valid_G) - rho_G * g_lag(valid_G);
sigma_G = std(eps_G, "omitnan");

%% 3. Collect results

ar_calibration = table( ...
    ["rho_A"; "sigma_A"; "rho_G"; "sigma_G"], ...
    [rho_A; sigma_A; rho_G; sigma_G], ...
    'VariableNames', {'Parameter', 'Value'} ...
);

disp(ar_calibration);
Output
Output
    Parameter     Value  
    _________    ________

    "rho_A"       0.68913
    "sigma_A"    0.012113
    "rho_G"       0.56783
    "sigma_G"    0.017086


The estimated technology process is moderately persistent, with $\rho_A \approx 0.69$. The government spending cycle is also persistent, but somewhat less so, with $\rho_G \approx 0.57$. The innovation standard deviations imply that the government spending shock is slightly more volatile than the technology shock in this calibration.

6. Steady-state calibration

The model is normalized with $A=1$ and $N=1/3$. Given the empirical ratios, the code computes the implied steady-state values of output, consumption, investment, government spending, capital, wages, tax rate, interest rate, discount factor and the preference weight on leisure.

Dynare model codeDynare
Code
var c invest k y n w r tau a g;
varexo eps_a eps_g;

parameters beta alpha delta theta rho_a rho_g sigma_a sigma_g g_ss;

beta    = 0.9718;
alpha   = 0.46718;
delta   = 0.036346;
theta   = 1.2786;

rho_a   = 0.68913;
sigma_a = 0.012113;

rho_g   = 0.56783;
sigma_g = 0.017086;

g_ss    = 0.30512;


model;

    // Euler equation
    1/c = beta * (1/c(+1)) * (1 + r(+1));

    // Labor supply
    w = theta/(1 - tau) * c/(1 - n);

    // Labor demand
    w = (1 - alpha) * y / n;

    // Rental rate of capital
    r = alpha * y / k(-1) - delta;

    // Government budget constraint
    tau = g / (w * n);

    // Capital accumulation
    k = (1 - delta) * k(-1) + invest;

    // Resource constraint
    y = c + invest + g;

    // Production function
    y = a * k(-1)^alpha * n^(1 - alpha);

    // Technology process
    log(a) = rho_a * log(a(-1)) + eps_a;

    // Government spending process
    log(g) = (1 - rho_g) * log(g_ss) + rho_g * log(g(-1)) + eps_g;

end;

initval;
    a      = 1;
    n      = 0.33333;
    y      = 1.8697;
    c      = 1.0810;
    invest = 0.48314;
    g      = 0.30512;
    k      = 13.363;
    w      = 2.9887;
    tau    = 0.30628;
    r      = 0.029022;
end;


resid;
steady;
check;


shocks;
    var eps_a; stderr 0.012113;
    var eps_g; stderr 0.017086;
end;


stoch_simul(order = 1, irf = 40, periods=20000, drop=1000)
    y c invest k n w r tau a g;
MATLAB codeCell 5
Code
%% =========================================================
%  Full RBC calibration table
%  Normalization: A = 1, N = 1/3
%  Y is implied by the production function
%% =========================================================

%% 1. Ratios from data

C_Y = italy_rbc.C ./ italy_rbc.Y;
I_Y = italy_rbc.I ./ italy_rbc.Y;
G_Y = italy_rbc.G ./ italy_rbc.Y;
K_Y = italy_rbc.K ./ italy_rbc.Y;

delta_series = italy_rbc.I ./ italy_rbc.K;
alpha_series = italy_rbc.alpha;

c_y_bar   = mean(C_Y, "omitnan");
i_y_bar   = mean(I_Y, "omitnan");
g_y_bar   = mean(G_Y, "omitnan");
k_y_bar   = mean(K_Y, "omitnan");
delta_bar = mean(delta_series, "omitnan");
alpha_bar = mean(alpha_series, "omitnan");


A_ss = 1;
N_ss = 1/3;

alpha_ss = alpha_bar;
delta_ss = delta_bar;

% From production function:
% Y = A * K^alpha * N^(1-alpha)
% K = (K/Y) * Y
%
% Therefore:
% Y^(1-alpha) = A * (K/Y)^alpha * N^(1-alpha)

Y_ss = (A_ss * k_y_bar^alpha_ss * N_ss^(1 - alpha_ss))^(1 / (1 - alpha_ss));

K_ss = k_y_bar * Y_ss;
C_ss = c_y_bar * Y_ss;
I_ss = i_y_bar * Y_ss;
G_ss = g_y_bar * Y_ss;


% Real wage
w_ss = (1 - alpha_ss) * Y_ss / N_ss;

% Labor income tax
tau_ss = G_ss / (w_ss * N_ss);
% equivalently:
% tau_ss = g_y_bar / (1 - alpha_ss);

% Real interest rate
r_ss = alpha_ss * Y_ss / K_ss - delta_ss;

% Discount factor
beta_ss = 1 / (1 + r_ss);

% Preference weight on leisure
theta_ss = w_ss * (1 - tau_ss) * (1 - N_ss) / C_ss;

if exist("rho_A", "var") ~= 1
    rho_A = NaN;
end

if exist("sigma_A", "var") ~= 1
    sigma_A = NaN;
end

if exist("rho_G", "var") ~= 1
    rho_G = NaN;
end

if exist("sigma_G", "var") ~= 1
    sigma_G = NaN;
end


full_calibration = table( ...
    ["A"; "N"; "Y"; "C"; "I"; "G"; "K"; ...
     "C/Y"; "I/Y"; "G/Y"; "K/Y"; ...
     "alpha"; "delta"; "w"; "tau"; "r"; "beta"; "theta"; ...
     "rho_A"; "sigma_A"; "rho_G"; "sigma_G"], ...
    [A_ss; N_ss; Y_ss; C_ss; I_ss; G_ss; K_ss; ...
     c_y_bar; i_y_bar; g_y_bar; k_y_bar; ...
     alpha_ss; delta_ss; w_ss; tau_ss; r_ss; beta_ss; theta_ss; ...
     rho_A; sigma_A; rho_G; sigma_G], ...
    'VariableNames', {'Object', 'Value'} ...
);

disp(full_calibration);

%% 6. Checks

production_gap = Y_ss - A_ss * K_ss^alpha_ss * N_ss^(1 - alpha_ss);
resource_gap   = Y_ss - C_ss - I_ss - G_ss;

fprintf("\nProduction function check: Y - A*K^alpha*N^(1-alpha) = %.10f\n", production_gap);
fprintf("Resource constraint check: Y - C - I - G = %.10f\n", resource_gap);

% If resource_gap is not zero, this is the implied net export / residual share
NX_Y = resource_gap / Y_ss;
fprintf("Implied residual share, (Y - C - I - G)/Y = %.6f\n", NX_Y);
Output
Output
     Object       Value  
    _________    ________

    "A"                 1
    "N"           0.33333
    "Y"            1.8697
    "C"             1.081
    "I"           0.48314
    "G"           0.30512
    "K"            13.363
    "C/Y"         0.57816
    "I/Y"          0.2584
    "G/Y"         0.16319
    "K/Y"          7.1469
    "alpha"       0.46718
    "delta"      0.036346
    "w"            2.9887
    "tau"         0.30628
    "r"          0.029022
    "beta"         0.9718
    "theta"        1.2786
    "rho_A"       0.68913
    "sigma_A"    0.012113
    "rho_G"       0.56783
    "sigma_G"    0.017086


Production function check: Y - A*K^alpha*N^(1-alpha) = 0.0000000000
Resource constraint check: Y - C - I - G = 0.0004717935
Implied residual share, (Y - C - I - G)/Y = 0.000252

The steady state is internally consistent: the production function check is essentially zero, and the resource constraint residual is very small. The implied labor income tax rate is about 0.31, which reflects the government spending share relative to labor income. The small residual share indicates that the calibrated closed-economy resource constraint is a good approximation for this exercise.

7. Business-cycle moments

The next block computes empirical and model-implied business-cycle statistics for the same variables: $Y$, $C$, $I$, $N$, $Y/N$, $A$ and $G$. The empirical series are log-transformed and HP-filtered. The model series are simulated with Dynare and treated in the same way, so the reported moments are comparable.

MATLAB codeCell 6
Code
cd('C:\Users\User\Documents\Dynamic_macro\RBC project')

clear
clc
close all

lambda_emp = 100;
lambda_model = 100;

load("italy_rbc_2010_2023.mat", "italy_rbc")

emp = struct();

emp.Y = italy_rbc.Y;
emp.C = italy_rbc.C;
emp.I = italy_rbc.I;
emp.N = italy_rbc.N;
emp.YN = italy_rbc.Y ./ italy_rbc.N;
emp.A = italy_rbc.A;
emp.G = italy_rbc.G;

emp_var_names = ["Y", "C", "I", "N", "YN", "A", "G"];
emp_row_names = ["Y", "C", "I", "N", "Y/N", "A", "G"];

nvars = numel(emp_var_names);
cycles = cell(nvars, 1);

for j = 1:nvars

    vname = emp_var_names(j);
    x = emp.(vname);
    x = x(:);

    z = log(x);

    cyc = NaN(size(z));
    valid = ~isnan(z) & isfinite(z);
    zv = z(valid);
    Tn = length(zv);

    if Tn >= 4
        Imat = speye(Tn);
        D = spdiags([ones(Tn-2,1), -2*ones(Tn-2,1), ones(Tn-2,1)], 0:2, Tn-2, Tn);
        trend = (Imat + lambda_emp * (D' * D)) \ zv;
        cyc(valid) = zv - trend;
    end

    cycles{j} = 100 * cyc;

end

y_idx = find(emp_var_names == "Y", 1);
y_cyc = cycles{y_idx};

std_dev = NaN(nvars, 1);
rel_std = NaN(nvars, 1);
ac1 = NaN(nvars, 1);
corr_y = NaN(nvars, 1);

std_y = std(y_cyc(~isnan(y_cyc) & isfinite(y_cyc)));

for j = 1:nvars

    x = cycles{j};

    valid_std = ~isnan(x) & isfinite(x);
    std_dev(j) = std(x(valid_std));
    rel_std(j) = std_dev(j) / std_y;

    x_now = x(2:end);
    x_lag = x(1:end-1);

    valid_ac = ~isnan(x_now) & ~isnan(x_lag) & isfinite(x_now) & isfinite(x_lag);

    if sum(valid_ac) > 2
        ac1(j) = corr(x_now(valid_ac), x_lag(valid_ac));
    end

    valid_corr = ~isnan(x) & ~isnan(y_cyc) & isfinite(x) & isfinite(y_cyc);

    if sum(valid_corr) > 2
        corr_y(j) = corr(x(valid_corr), y_cyc(valid_corr));
    end

end

empirical_table = table( ...
    emp_row_names(:), ...
    std_dev, ...
    rel_std, ...
    corr_y, ...
    'VariableNames', { ...
        'Variable', ...
        'StandardDeviation', ...
        'RelativeStandardDeviation', ...
        'ContemporaneousCorrelationWithOutput' ...
    } ...
);

set(0, 'DefaultFigureVisible', 'off')
dynare rbc_project.mod noclearall
close all
set(0, 'DefaultFigureVisible', 'on')

if ischar(M_.endo_names)
    endo_names = string(cellstr(M_.endo_names));
elseif iscell(M_.endo_names)
    endo_names = string(M_.endo_names);
else
    endo_names = string(M_.endo_names);
end

endo_names = strtrim(endo_names(:));

idx_y = find(endo_names == "y", 1);
idx_c = find(endo_names == "c", 1);
idx_i = find(endo_names == "invest", 1);
idx_n = find(endo_names == "n", 1);
idx_a = find(endo_names == "a", 1);
idx_g = find(endo_names == "g", 1);

idx_all = [idx_y idx_c idx_i idx_n idx_a idx_g];

if any(isempty(idx_all)) || numel(idx_all) < 6
    disp(endo_names)
    error("Hianyzik legalabb egy valtozo a Dynare modellbol.")
end

sim_raw = oo_.endo_simul;
idx_max = max(idx_all);

if size(sim_raw, 1) >= idx_max && size(sim_raw, 2) > 10
    y_sim = sim_raw(idx_y, :)';
    c_sim = sim_raw(idx_c, :)';
    i_sim = sim_raw(idx_i, :)';
    n_sim = sim_raw(idx_n, :)';
    a_sim = sim_raw(idx_a, :)';
    g_sim = sim_raw(idx_g, :)';
elseif size(sim_raw, 2) >= idx_max && size(sim_raw, 1) > 10
    y_sim = sim_raw(:, idx_y);
    c_sim = sim_raw(:, idx_c);
    i_sim = sim_raw(:, idx_i);
    n_sim = sim_raw(:, idx_n);
    a_sim = sim_raw(:, idx_a);
    g_sim = sim_raw(:, idx_g);
else
    disp(size(sim_raw))
    error("Nincs eleg szimulalt megfigyeles. A .mod fajlban legyen periods = 20000.")
end

model = struct();

model.Y = y_sim;
model.C = c_sim;
model.I = i_sim;
model.N = n_sim;
model.YN = y_sim ./ n_sim;
model.A = a_sim;
model.G = g_sim;

model_var_names = ["Y", "C", "I", "N", "YN", "A", "G"];
model_row_names = ["Y", "C", "I", "N", "Y/N", "A", "G"];

nvars = numel(model_var_names);
cycles = cell(nvars, 1);

for j = 1:nvars

    vname = model_var_names(j);
    x = model.(vname);
    x = x(:);

    z = log(x);

    cyc = NaN(size(z));
    valid = ~isnan(z) & isfinite(z);
    zv = z(valid);
    Tn = length(zv);

    if Tn >= 4
        Imat = speye(Tn);
        D = spdiags([ones(Tn-2,1), -2*ones(Tn-2,1), ones(Tn-2,1)], 0:2, Tn-2, Tn);
        trend = (Imat + lambda_model * (D' * D)) \ zv;
        cyc(valid) = zv - trend;
    end

    cycles{j} = 100 * cyc;

end

y_idx = find(model_var_names == "Y", 1);
y_cyc = cycles{y_idx};

std_dev = NaN(nvars, 1);
rel_std = NaN(nvars, 1);
ac1 = NaN(nvars, 1);
corr_y = NaN(nvars, 1);

std_y = std(y_cyc(~isnan(y_cyc) & isfinite(y_cyc)));

for j = 1:nvars

    x = cycles{j};

    valid_std = ~isnan(x) & isfinite(x);
    std_dev(j) = std(x(valid_std));
    rel_std(j) = std_dev(j) / std_y;

    x_now = x(2:end);
    x_lag = x(1:end-1);

    valid_ac = ~isnan(x_now) & ~isnan(x_lag) & isfinite(x_now) & isfinite(x_lag);

    if sum(valid_ac) > 2
        ac1(j) = corr(x_now(valid_ac), x_lag(valid_ac));
    end

    valid_corr = ~isnan(x) & ~isnan(y_cyc) & isfinite(x) & isfinite(y_cyc);

    if sum(valid_corr) > 2
        corr_y(j) = corr(x(valid_corr), y_cyc(valid_corr));
    end

end

model_table = table( ...
    model_row_names(:), ...
    std_dev, ...
    rel_std, ...
    corr_y, ...
    'VariableNames', { ...
        'Variable', ...
        'StandardDeviation', ...
        'RelativeStandardDeviation', ...
        'ContemporaneousCorrelationWithOutput' ...
    } ...
);
Output
Output
Starting Dynare (version 6.5).
Calling Dynare with arguments: noclearall
Starting preprocessing of the model file ... 
Found 10 equation(s). 
Evaluating expressions... 
Computing static model derivatives (order 1). 
Normalizing the static model... 
Finding the optimal block decomposition of the static model... 
3 block(s) found: 
  2 recursive block(s) and 1 simultaneous block(s). 
  the largest simultaneous block has 8 equation(s) 
                                 and 8 feedback variable(s). 
Computing dynamic model derivatives (order 1). 
Normalizing the dynamic model... 
Finding the optimal block decomposition of the dynamic model... 
3 block(s) found: 
  2 recursive block(s) and 1 simultaneous block(s). 
  the largest simultaneous block has 8 equation(s) 
                                 and 7 feedback variable(s). 
Preprocessing completed. 
Preprocessing time: 0h00m01s.


Residuals of the static equations:
Equation number  1: 1   :      -0.000003
Equation number  2: w   :       0.000117
Equation number  3: 3   :       0.000029
Equation number  4: r   :       0.000002
Equation number  5: tau :       0.000003
Equation number  6: k   :       0.002552
Equation number  7: y   :       0.000440
Equation number  8: 8   :      -0.000025
Equation number  9: 9   :       0.000000
Equation number 10: 10  :       0.000000



STEADY-STATE RESULTS:

c      		 1.08101
invest 		 0.486455
k      		 13.384
y      		 1.87259
n      		 0.333829
w      		 2.98882
r      		 0.0290183
tau    		 0.305807
a      		 1
g      		 0.30512

EIGENVALUES:
         Modulus             Real        Imaginary
          0.5678           0.5678                0
          0.6891           0.6891                0
          0.9527           0.9527                0
           1.101            1.101                0
       1.698e+18       -1.698e+18                0

There are 2 eigenvalue(s) larger than 1 in modulus for 2 forward-looking variable(s)
The rank condition is verified.


MODEL SUMMARY

  Number of variables:         10
  Number of stochastic shocks: 2
  Number of state variables:   3
  Number of jumpers:           2
  Number of static variables:  5


MATRIX OF COVARIANCE OF EXOGENOUS SHOCKS
Variables       eps_a     eps_g
eps_a        0.000147  0.000000
eps_g        0.000000  0.000292

POLICY AND TRANSITION FUNCTIONS
                                   y               c          invest               k               n               w               r             tau               a               g
Constant                    1.872588        1.081013        0.486455       13.383991        0.333829        2.988815        0.029018        0.305807        1.000000        0.305120
k(-1)                       0.053121        0.064087       -0.010965        0.952689       -0.004096        0.121461       -0.003030       -0.008675               0               0
a(-1)                       2.316885        0.257332        2.059553        2.059553        0.343423        0.623236        0.080873       -0.378364        0.689130               0
g(-1)                      -0.912094       -0.161304       -1.318620       -1.318620       -0.305169        1.276439       -0.031837        0.718061               0        0.567830
eps_a                       3.362043        0.373415        2.988627        2.988627        0.498343        0.904381        0.117355       -0.549046        1.000000               0
eps_g                      -0.490108       -0.086676       -0.708552       -0.708552       -0.163981        0.685887       -0.017108        0.385846               0        0.305120


MOMENTS OF SIMULATED VARIABLES
VARIABLE              MEAN       STD. DEV.        VARIANCE        SKEWNESS        KURTOSIS
y                 1.870602        0.064902        0.004212        0.021178        0.026219
c                 1.079676        0.026287        0.000691        0.017035        0.197816
invest            0.485753        0.050572        0.002558        0.014586       -0.010129
k                13.364886        0.370965        0.137615        0.017047        0.206034
n                 0.333747        0.008611        0.000074        0.011088       -0.009628
w                 2.986376        0.052745        0.002782        0.013068        0.158308
r                 0.029042        0.001981        0.000004        0.002137       -0.009371
tau               0.306184        0.013117        0.000172       -0.001072        0.029705
a                 0.999736        0.016577        0.000275        0.022028       -0.019115
g                 0.305173        0.006195        0.000038       -0.033179       -0.009730


CORRELATION OF SIMULATED VARIABLES
VARIABLE          y        c   invest        k        n        w        r      tau        a        g
y            1.0000   0.7273   0.9250   0.6532   0.8698   0.6926   0.6433  -0.8842   0.9424  -0.1608
c            0.7273   1.0000   0.4255   0.9947   0.3322   0.9428  -0.0571  -0.6336   0.4860  -0.0970
invest       0.9250   0.4255   1.0000   0.3299   0.9898   0.3699   0.8740  -0.8792   0.9564  -0.2785
k            0.6532   0.9947   0.3299   1.0000   0.2339   0.9409  -0.1594  -0.5613   0.3976  -0.0708
n            0.8698   0.3322   0.9898   0.2339   1.0000   0.2465   0.9007  -0.8815   0.9092  -0.3774
w            0.6926   0.9428   0.3699   0.9409   0.2465   1.0000  -0.0530  -0.4481   0.5220   0.2358
r            0.6433  -0.0571   0.8740  -0.1594   0.9007  -0.0530   1.0000  -0.5921   0.8241  -0.1527
tau         -0.8842  -0.6336  -0.8792  -0.5613  -0.8815  -0.4481  -0.5921   1.0000  -0.7599   0.6033
a            0.9424   0.4860   0.9564   0.3976   0.9092   0.5220   0.8241  -0.7599   1.0000   0.0033
g           -0.1608  -0.0970  -0.2785  -0.0708  -0.3774   0.2358  -0.1527   0.6033   0.0033   1.0000


AUTOCORRELATION OF SIMULATED VARIABLES
VARIABLE         1       2       3       4       5
y           0.7679  0.6047  0.4973  0.4200  0.3617
c           0.9823  0.9565  0.9262  0.8929  0.8579
invest      0.6603  0.4316  0.2883  0.1918  0.1233
k           0.9910  0.9707  0.9433  0.9117  0.8777
n           0.6389  0.4004  0.2530  0.1557  0.0870
w           0.9424  0.8941  0.8524  0.8139  0.7776
r           0.6579  0.4245  0.2777  0.1780  0.1078
tau         0.7024  0.5146  0.3993  0.3245  0.2702
a           0.6845  0.4651  0.3255  0.2289  0.1610
g           0.5593  0.3061  0.1652  0.0873  0.0429



VARIANCE DECOMPOSITION SIMULATING ONE SHOCK AT A TIME (in percent)
                     eps_a             eps_g  Tot. lin. contr.
y                    96.29              3.22             99.52
c                    93.10              5.14             98.24
invest               92.36              7.93            100.28
k                    93.10              5.10             98.19
n                    85.87             14.52            100.39
w                    89.02              8.74             97.76
r                    96.55              3.85            100.41
tau                  62.87             37.25            100.12
a                   100.01              0.00            100.01
g                     0.00            100.01            100.01
Note: numbers do not add up to 100 due to non-zero correlation of simulated shocks in small samples

Total computing time : 0h00m06s

8. Model versus data

The final block prints the two tables: first the simulated RBC moments, then the empirical moments from the Italian data. The comparison focuses on standard deviations, relative standard deviations and contemporaneous correlations with output.

MATLAB codeCell 7
Code
fprintf("\nModel business-cycle statistics\n")
disp(model_table)

fprintf("\nEmpirical business-cycle statistics\n")
disp(empirical_table)
Output
Output
Model business-cycle statistics
    Variable    StandardDeviation    RelativeStandardDeviation    ContemporaneousCorrelationWithOutput
    ________    _________________    _________________________    ____________________________________

     "Y"              2.0342                        1                                  1              
     "C"             0.45503                  0.22368                            0.75725              
     "I"              7.4142                   3.6447                            0.98661              
     "N"              1.8621                  0.91535                            0.96715              
     "Y/N"           0.52774                  0.25943                             0.4422              
     "A"              1.1327                  0.55683                            0.97754              
     "G"               1.566                  0.76981                           -0.17492              


Empirical business-cycle statistics
    Variable    StandardDeviation    RelativeStandardDeviation    ContemporaneousCorrelationWithOutput
    ________    _________________    _________________________    ____________________________________

     "Y"              3.2616                        1                                  1              
     "C"              3.3666                   1.0322                            0.91908              
     "I"              8.3301                    2.554                            0.87074              
     "N"              3.5901                   1.1007                            0.86986              
     "Y/N"            1.7764                  0.54466                           0.078091              
     "A"             0.90964                  0.27889                            0.82819              
     "G"              2.2525                   0.6906                            0.49963              


The model reproduces some qualitative RBC patterns, especially the high procyclicality of investment, labor and productivity. Investment is much more volatile than output in the model, and it is also highly correlated with output, which is consistent with the empirical table.

There are also clear mismatches. Model consumption is far too smooth relative to output: its relative standard deviation is about 0.22 in the model, compared with about 1.03 in the data. The model also makes productivity much more tightly correlated with output than in the data. Government spending is the largest discrepancy: empirically it is positively correlated with output, while in the model it is slightly negatively correlated. This happens because government spending is imposed as an exogenous process rather than modeled as a systematic fiscal response to the cycle.

Overall, the calibrated RBC model captures the basic technology-driven comovement mechanism, but it does not fully match the volatility and fiscal comovement patterns of the Italian data.