% This program replicates tables A.2, A.3, and A.4 in Online Appendix of
% "A Static Capital Buffer is Hard to Beat" by Matthew Canzoneri,
% Behzad Diba, Luca Guerrieri and Arsenii Mishin.
%
% This program computes the results (parameter values of the shock processes,
% variance decomposition, and matching moments) using a SMM procedure under
% an alternative calibration that considers TFP and ISP shocks (there
% is no volatility shock compared with the main framework)


clear all

setpath_windows

global M_ oo_ options_
global overwrite overwrite_param_names

%% set some switches

short_sample_switch = 0; % use a sample from 1980 through 2008; if set to zero, use the benchmark sample from 1980 to 2024
plot_irf_switch = 0;     % switch to plot IRFs
overwrite_results = 0;   % set to 1 to overwrite the results used for the paper
load_from_disk = 1;      % set to 1 to load initial guess from disk

%% model setup
mod00 = 'model_safe_all_shocks';
mod11 = 'model_mixed_all_shocks';
mod10 = 'model_risky_all_shocks';
mod01 = 'model_spare_all_shocks';

constraint1 = 'chi2_upper_bar+chi2_upper_bar_ss<0';    % If the Lagrange multiplier on the
% no-shorting constraints for risky
% loans is slack, it means that I want
% to make risky loans.
% Accordingly, go into 11 or 10

constraint_relax1 = '((l+l_ss)<0)|(l_upper_bar<0)'; % If it turns out that the quantity of
% risky loans is negative,  00
% notice that l is set to 0 in the 10
% file, therefore the first condition
% is irrelevant when transition back
% from 10, but it might be relevant
% when transition back from 11.


constraint2 = '((chi2_lower_bar<0)&(l_upper_bar>0))'; % If the Lagrange multiplier on the
% no-shorting constraints for risky
% loans  and on safe loans are both slack, it means that I want
% to make risky and safe loans. Go to
% 01 or 11, but since if this
% condition is met constraint1 will
% also be active, we can only go into
% 11. Hence, 01 is spare.

constraint_relax2 ='((l+l_ss)<0)|(l_upper_bar<0)';  % return to 10 or 00, but
% notice that
% constraint_relax1 is identical,
% therefore, if either
% condition is violated
% we will go back to 00


%% This part does some setup and error checking
% We want to make sure that two functions yield the same results
% solve_two_constraints_add_violvecbool
% solve_two_constraints_add_violvecbool_noparse
% the _noparse function will be used in the moment matching routines --- it
% excludes error checking on the .mod files and does not parse the .mod files. That's why
% solve_two_constraints_violvecbool needs to be run first.

nperiods = 100;
maxiter = 20;




% select shock for IRFs
% irfshock = char('eps_gamma');
% shockssequence = -0.01;

irfshock = char('eps_a');
shockssequence = randn(nperiods,1);

% % load monthly data for consumption excluding food and energy, turn it to
% % quarterly by averaging
% data_table_m = readtable('Commercial_banks.xlsx','Sheet',"monthly");
% data_table_m = table2timetable(data_table_m);
% data_table_q = retime(data_table_m,'quarterly','mean');
% data_table_q = timetable2table(data_table_q);
% 
% data_table_q = data_table_q(year(data_table_q.observation_date)>1979 &...
%     year(data_table_q.observation_date)<2025,"DPCCRA3M086SBEA");

% import other quarterly data
data_table = readtable('Commercial_banks.xlsx','Sheet',"Quarterly");
data_table = data_table(year(data_table.observation_date)>1979 & year(data_table.observation_date)<2025,:);


%dow_return = log(data_table.dow(2:end))-log(data_table.dow(1:end-1));
%dow_return_detrended = 100*bpass(dow_return,6,32);

% detrend data
gdp_detrended = 100*bpass(log(data_table.A191RA3Q086SBEA),6,32);

investment_detrended = 100*bpass(log(data_table.A007RA3Q086SBEA),6,32);

consumption_detrended = 100*bpass(log(data_table.DPCCRA3Q086SBEA),6,32);

investment_price_detrended = 100*bpass(log(data_table.B007RG3Q086SBEA)-log(data_table.JCXFE),6,32);

% assemble series for moment matching
% dataset = [gdp_detrended(2:end) investment_detrended(2:end) consumption_detrended(2:end)];
% var_labels =char('GDP','Investment','Consumption');

dataset = [gdp_detrended(2:end) investment_detrended(2:end) investment_price_detrended(2:end)];
var_labels =char('GDP','Investment','Investment Price');

% select strategy for creating the GMM weighting matrix
%optwt =1; % weighting matrix is diagonal optimal SMM matrix
optwt =2; % optimal GMM matrix
%optwt = -1; % weighting matrix is the identity matrix

% create the GMM weighting matrix
[data_momm,wtmat] = wtm_rm_noskew_opt(dataset,optwt);

%'sigma_eps_taup'
%overwrite_param_names = char('taup');
%overwrite = [0.025];
%overwrite_param_names = char('xip');
%overwrite = 0.003;


% solve the model
% eval(['dynare ',modnam,' noclearall']
[ zdatalinear_, zdatapiecewise1_, zdatass_, oo00_,  M00_, violvecbool,...
    M10_,M01_,M11_,oo10_,oo01_,oo11_,...
    options00_,options10_,options01_,options11_] = ...
    solve_two_constraints_add_violvecbool(mod00,mod10,mod01,mod11,...
    constraint1, constraint2,...
    constraint_relax1, constraint_relax2,...
    shockssequence,irfshock,nperiods,0,maxiter);

[zdatapiecewise2_, zdatass_, oo00_,  M00_, violvecbool_last_shock] = ...
    solve_two_constraints_add_violvecbool_no_parse(M00_,M10_,M01_,M11_,oo00_,...
    options00_,constraint1, constraint2,...
    constraint_relax1, constraint_relax2,...
    shockssequence,irfshock,nperiods,0,maxiter);

nperiods = 50;
[zdatapiecewise3_, zdatass_, oo00_,  M00_, violvecbool_last_shock] = ...
    solve_two_constraints_add_violvecbool_no_parse_short(M00_,M10_,M01_,M11_,oo00_,...
    options00_,constraint1, constraint2,...
    constraint_relax1, constraint_relax2,...
    shockssequence,irfshock,nperiods,0,maxiter);



if max(max(abs(zdatapiecewise1_-zdatapiecewise2_)))>10^(-15)
    error('Check no_parse')
end

%% set plot_irf_switch=1 at the beginning of the program to plot IRFs
if plot_irf_switch

    prefix1 = 's1_';
    prefix2 = 's2_';
    for this_var = 1:M00_.endo_nbr
        eval([prefix1,deblank(M00_.endo_names{this_var,:}),'_irf=zdatapiecewise1_(:,this_var);'])
        eval([prefix1,deblank(M00_.endo_names{this_var,:}),'_ss=zdatass_(this_var);'])

        eval([prefix2,deblank(M00_.endo_names{this_var,:}),'_irf=zdatapiecewise3_(:,this_var);'])
        eval([prefix2,deblank(M00_.endo_names{this_var,:}),'_ss=zdatass_(this_var);'])

    end

    spread_title = 'Expected Equity Return Spread';
    spread_title = [spread_title newline '(risky-safe)'];

    titlelist = char('Bank capital requirement','Total deposits',...
        'Consumption', 'Total output',...
        spread_title, 'Total equity in banks' ,...
        'Expected safe equity return','Total capital' );


    line1 = 100*[ s1_gamma_var_irf, (s1_d_irf+s1_d_upper_bar_irf)/(s1_d_ss+s1_d_upper_bar_ss),...
        s1_c_irf/s1_c_ss,  s1_ytot_irf/s1_ytot_ss,...
        [s1_r_e_upper_bar_irf(2:end)+s1_r_e_upper_bar_ss-s1_r_e_irf(2:end)-s1_r_e_ss; nan], (s1_e_irf+s1_e_upper_bar_irf)/(s1_e_ss+s1_e_upper_bar_ss),...
        [s1_r_e_irf(2:end); nan], s1_ktot_irf/s1_ktot_ss ];

    line2 = 100*[ s2_gamma_var_irf, (s2_d_irf+s2_d_upper_bar_irf)/(s2_d_ss+s2_d_upper_bar_ss),...
        s2_c_irf/s2_c_ss,  s2_ytot_irf/s2_ytot_ss,...
        [s2_r_e_upper_bar_irf(2:end)+s2_r_e_upper_bar_ss-s2_r_e_irf(2:end)-s2_r_e_ss; nan], (s2_e_irf+s2_e_upper_bar_irf)/(s2_e_ss+s2_e_upper_bar_ss),...
        [s2_r_e_irf(2:end); nan], s2_ktot_irf/s2_ktot_ss ];


    percent = 'Percent';
    ppt = 'Perc. Point';
    level = 'Level';
    ylabels = char(ppt,percent,...
        percent,percent,...
        ppt,percent,...
        ppt,percent);
    figlabel = '';



    legendlist  = char('Higher capital requirements','Lower capital requirements');
    makechart(titlelist,legendlist,figlabel,ylabels,line1,line2);
end


%% moment matching

if load_from_disk
    % load starting guess from disk
    load my_estimates_constraints_spreads_my;
    init_theta = distance_param;
else
    %specify starting guess (only necessary the first time the codes are
    %run)
init_theta = [0.010112136254779
   0.009788591910459
   0.001651443864956
   0.955918339672063
   0.632636489829951
   0.709401893328523
   0.173765894644399
   0.761639078395671
   0.001207557515828
   0.046257298913245];


init_theta = [0.010112136254779
   0.009788591910459
   0.001651443864956
   0.955918339672063
   0.632636489829951
   0.709401893328523
   0.173765894644399
   0.761639078395671
   0.001207557515828
   0.046257298913243];
% 



end

% parameters to be calibrated
overwrite_param_names = char('sigma_eps_a','sigma_eps_isp','sigma_eps_tau',...
    'rhop_a','rho_isp','rho_taup','phip_i','kappap','xip','taup');

% test the distance function
[dist, model_momm, average_failure_rate, average_gdp_loss, mean_negative_equity_at_default] = distance_function_constraints_spreads(init_theta,data_momm,wtmat,...
    M00_,M10_,M01_,M11_,oo00_,...
    options00_,constraint1, constraint2,...
    constraint_relax1, constraint_relax2);

mydist = @(init_theta) distance_function_constraints_spreads(init_theta,data_momm,wtmat,...
M00_,M10_,M01_,M11_,oo00_,...
options00_,constraint1, constraint2,...
constraint_relax1, constraint_relax2);


% options for the numerical minimizer
optimization_options = optimset('display','iter','MaxIter',150,'MaxFunEvals',1e10,'TolFun',10e-3,'TolX',1e-4);


lowerbound = [0.0001
              0.0001
              0.0001
              0
              0
              0
              0
              0
              0
              0];

upperbound = [1
    1
    1
    .999
    .999
    .999
    200
    .95
    1
    .1];

%thetaunc=simulannealbnd(mydist,init_theta,lowerbound,upperbound,optimization_options);
%thetaunc=particleswarm(mydist,init_theta,lowerbound,upperbound,optimization_options);

options_patternsearch = optimoptions('patternsearch');
options_patternsearch = optimoptions(options_patternsearch,'display','iter','MaxIterations',100,'MeshTolerance',1e-5,'StepTolerance',1e-5);

% iterate through guesses -- set a number for the rep loop; fminsearch has
% better convergence properties if we set a low maxfunevals (above) and
% call it many times --- we have reset the loop to 1 iteration after
% ascertaining convergence
thetaunc = init_theta;
if overwrite_results
    for rep = 1:1
        rep

        thetaunc = patternsearch(mydist,init_theta,[],[],[],[],lowerbound,upperbound,[],options_patternsearch);
        init_theta=thetaunc;
        distance_param = thetaunc;
        save my_estimates_constraints_spreads_my distance_param rep
        
        % [thetaunc,fval,exitflag] = ...
        %     fminsearch(mydist,init_theta,optimization_options);
        % init_theta = thetaunc;
        % distance_param = thetaunc;
        % save my_estimates_constraints_spreads_my distance_param rep

       
    end

    thetaunc = patternsearch(mydist,init_theta,[],[],[],[],lowerbound,upperbound,[],options_patternsearch);
        init_theta=thetaunc;
        distance_param = thetaunc;
        save my_estimates_constraints_spreads_my distance_param rep

    [dist, model_momm, average_failure_rate, average_gdp_loss] = distance_function_constraints_spreads(init_theta,data_momm,wtmat,...
    M00_,M10_,M01_,M11_,oo00_,...
    options00_,constraint1, constraint2,...
    constraint_relax1, constraint_relax2);
end


%% display moments
if isnan(model_momm)~=1
    [model_momm data_momm]

    nvars = size(var_labels,1);
    this_mom = 0;
    moment_labels = char('')
    for i = 1:nvars
        for j=i:nvars
            this_mom = this_mom+1;
            if i == j
                moment_labels = strvcat(moment_labels,['Var(',deblank(var_labels(i,:)),')']);
                display(['Var(',deblank(var_labels(i,:)),') in the model is ',num2str(model_momm(this_mom)),' and in the data is ',num2str(data_momm(this_mom))])
            else
                moment_labels = strvcat(moment_labels,['Corr(',deblank(var_labels(i,:)),',', deblank(var_labels(j,:)),')']);
                display(['The covariance of ',deblank(var_labels(i,:)),' with ',deblank(var_labels(j,:)),' in the model is ',num2str(model_momm(this_mom)),' and in the data is ',num2str(data_momm(this_mom))])

            end
        end
    end

    for i = 1:nvars
        this_mom = this_mom+1;
        moment_labels = strvcat(moment_labels,['Autocorr(',deblank(var_labels(i,:)),')']);
        display(['The auto-covariance of ',deblank(var_labels(i,:)),' is ',num2str(model_momm(this_mom)),' and in the data is ',num2str(data_momm(this_mom))])

    end

    display(['The average failure rate is ',num2str(average_failure_rate)])
    display(['The average GDP loss is ',num2str(average_gdp_loss)])
    display(['The mean negative equity at default is ',num2str(mean_negative_equity_at_default)])


    %% Shock decompositions for the different moments
    overwrite = init_theta;


    n_simul = 1000;

    irf_shock = char('eps_a','eps_isp','eps_tau');
    nshocks = size(irf_shock,1);
    nperiods = size(dataset,1);

    randn('seed',1);

    model_momm_mat = nan(size(model_momm,1),n_simul);
    mean_negative_equity_at_default_array = nan(n_simul,1);
    average_failure_rate_array = nan(n_simul,1);
    average_gdp_loss_array = nan(n_simul,1);
    max_failure_rate_array = nan(n_simul,1);
    


    for j=1:n_simul
        j
        shockssequence = randn(nperiods,nshocks);


        [zdatapiecewise, zdatass, oo00_,  M00_, violvecbool_last_shock] = ...
            solve_two_constraints_add_violvecbool_no_parse(M00_,M10_,M01_,M11_,oo00_,...
            options00_,constraint1, constraint2,...
            constraint_relax1, constraint_relax2,...
            shockssequence,irf_shock,nperiods,0,maxiter);


        [model_momm_mat(:,j),...
            mean_negative_equity_at_default_array(j),...
            average_failure_rate_array(j),...
            average_gdp_loss_array(j),...
            max_failure_rate_array(j),...
            crisis_prob(j)]=make_moments(M00_,zdatapiecewise,zdatass);

    end

    model_momm_mat= sort(model_momm_mat,2);
    mean_negative_equity_at_default_array = sort(mean_negative_equity_at_default_array);
    average_failure_rate_array = sort(average_failure_rate_array);
    average_gdp_loss_array = sort(average_gdp_loss_array);
    max_failure_rate_array = sort(max_failure_rate_array);
    crisis_prob = sort(crisis_prob);

    model_momm_lb = model_momm_mat(:,floor(n_simul*.05));
    model_momm_ub = model_momm_mat(:,ceil(n_simul*.95));
    average_failure_rate_lb = average_failure_rate_array(floor(n_simul*.05));
    average_failure_rate_ub = average_failure_rate_array(ceil(n_simul*.95));
    average_gdp_loss_lb = average_gdp_loss_array(floor(n_simul*.05));
    average_gdp_loss_ub = average_gdp_loss_array(ceil(n_simul*.95));



    sideways_switch = 0;
    caption = '';
    
    % hardwired first moments are computed in failure_rates_quarterly
    texfilename = 'table_moment';
    columnlabels = char('Data','Model','Model 5th perc.','Model 95th perc.');
    %rowlabels = moment_labels;
    rowlabels = char(moment_labels,'Avg. bank default rate, \%','Avg. GDP Gap with high default, \%');
    table_mat = [ data_momm                model_momm                model_momm_lb                model_momm_ub
                  400*(0.003799874586704)  400*average_failure_rate  400*average_failure_rate_lb  400*average_failure_rate_ub
                  -0.331394920309664       average_gdp_loss          average_gdp_loss_lb          average_gdp_loss_ub];

    column_precision = [2, 2, 2, 2];
    maketable(table_mat,texfilename,columnlabels,rowlabels,sideways_switch,caption,column_precision)




    nperiods = 5000;
    irf_shock = char('eps_a','eps_isp','eps_tau');

    randn('seed',1);
    nshocks = size(irf_shock,1);
    shockssequence = randn(nperiods,nshocks);


    [zdatapiecewise2_, zdatass_, oo00_,  M00_, violvecbool_last_shock] = ...
        solve_two_constraints_add_violvecbool_no_parse(M00_,M10_,M01_,M11_,oo00_,...
        options00_,constraint1, constraint2,...
        constraint_relax1, constraint_relax2,...
        shockssequence,irf_shock,nperiods,0,maxiter);

    ytot_pos = strmatch('ytot',M00_.endo_names,'exact');
    ltot2ytot_pos = strmatch('ltot2ytot',M00_.endo_names,'exact');
    gamma_var_pos = strmatch('gamma_var',M00_.endo_names,'exact');
    ginvest_pos = strmatch('ginvest',M00_.endo_names,'exact');
    q_pos = strmatch('q',M00_.endo_names,'exact');
    r_e_pos = strmatch('r_e',M00_.endo_names,'exact');

    f0_ytot_data = zdatapiecewise2_(:,ytot_pos);
    f0_ytot_ss = zdatass_(ytot_pos);
    f0_ltot2ytot_data = zdatapiecewise2_(:,ltot2ytot_pos);
    f0_ginvest_data = zdatapiecewise2_(:,ginvest_pos);
    f0_ginvest_ss = zdatass_(ginvest_pos);
    %f0_gamma_var_data = zdatapiecewise2_(:,gamma_var_pos);
    f0_q_data = zdatapiecewise2_(:,q_pos);
    f0_r_e_data = zdatapiecewise2_(:,r_e_pos);


    gdp = 100*bpass(f0_ytot_data/f0_ytot_ss,6,32);
    ltot2ytot = 100*bpass(f0_ltot2ytot_data,6,32);
    %gamma_var = 100*bpass(f0_gamma_var_data,6,32);
    ginvest = 100*bpass(f0_ginvest_data/f0_ginvest_ss,6,32);
    q = 100*bpass(f0_q_data,6,32);
    r_e =100*bpass(f0_r_e_data,6,32);

    tot_gdp_var = var(gdp);
    tot_ltot2ytot = var(ltot2ytot);
    %tot_gamma_var = var(gamma_var);
    tot_ginvest_var = var(ginvest);
    tot_q_var = var(q);
    tot_r_e_var = var(r_e);

    % repeat the variance calculation only turning on the TFP shock
    irf_shock = char('eps_a');

    [zdatapiecewise2_, zdatass_, oo00_,  M00_, violvecbool_last_shock] = ...
        solve_two_constraints_add_violvecbool_no_parse(M00_,M10_,M01_,M11_,oo00_,...
        options00_,constraint1, constraint2,...
        constraint_relax1, constraint_relax2,...
        shockssequence(:,1),irf_shock,nperiods,0,maxiter);

    ytot_pos = strmatch('ytot',M00_.endo_names,'exact');
    gamma_var_pos = strmatch('gamma_var',M00_.endo_names,'exact');


    f0_ytot_data = zdatapiecewise2_(:,ytot_pos);
    f0_ytot_ss = zdatass_(ytot_pos);
    f0_ltot2ytot_data = zdatapiecewise2_(:,ltot2ytot_pos);
    f0_ginvest_data = zdatapiecewise2_(:,ginvest_pos);
    f0_ginvest_ss = zdatass_(ginvest_pos);
    %f0_gamma_var_data = zdatapiecewise2_(:,gamma_var_pos);
    f0_q_data = zdatapiecewise2_(:,q_pos);
    f0_r_e_data = zdatapiecewise2_(:,r_e_pos);


    gdp = 100*bpass(f0_ytot_data/f0_ytot_ss,6,32);
    ltot2ytot = 100*bpass(f0_ltot2ytot_data,6,32);
    %gamma_var = 100*bpass(f0_gamma_var_data,6,32);
    ginvest = 100*bpass(f0_ginvest_data/f0_ginvest_ss,6,32);
    q = 100*bpass(f0_q_data,6,32);
    r_e =100*bpass(f0_r_e_data,6,32);

    part_gdp_var_a = var(gdp);
    part_ltot2ytot_var_a = var(ltot2ytot);
    %part_gamma_var_a = var(gamma_var);
    part_ginvest_var_a = var(ginvest);
    part_q_var_a = var(q);
    part_r_e_var_a = var(r_e);




    % get simulated data with only the ISP shock turned on

    irf_shock = char('eps_isp');

    [zdatapiecewise2_, zdatass_, oo00_,  M00_, violvecbool_last_shock] = ...
        solve_two_constraints_add_violvecbool_no_parse(M00_,M10_,M01_,M11_,oo00_,...
        options00_,constraint1, constraint2,...
        constraint_relax1, constraint_relax2,...
        shockssequence(:,2),irf_shock,nperiods,0,maxiter);

    ytot_pos = strmatch('ytot',M00_.endo_names,'exact');
    gamma_var_pos = strmatch('gamma_var',M00_.endo_names,'exact');


    f0_ytot_data = zdatapiecewise2_(:,ytot_pos);
    f0_ytot_ss = zdatass_(ytot_pos);
    f0_ltot2ytot_data = zdatapiecewise2_(:,ltot2ytot_pos);
    f0_ginvest_data = zdatapiecewise2_(:,ginvest_pos);
    f0_ginvest_ss = zdatass_(ginvest_pos);
    %f0_gamma_var_data = zdatapiecewise2_(:,gamma_var_pos);
    f0_q_data = zdatapiecewise2_(:,q_pos);
    f0_r_e_data = zdatapiecewise2_(:,r_e_pos);


    gdp = 100*bpass(f0_ytot_data/f0_ytot_ss,6,32);
    ltot2ytot = 100*bpass(f0_ltot2ytot_data,6,32);
    %gamma_var = 100*bpass(f0_gamma_var_data,6,32);
    ginvest = 100*bpass(f0_ginvest_data/f0_ginvest_ss,6,32);
    q = 100*bpass(f0_q_data,6,32);
    r_e =100*bpass(f0_r_e_data,6,32);


    part_gdp_var_isp = var(gdp);
    part_ltot2ytot_var_isp = var(ltot2ytot);
    %part_gamma_var_isp = var(gamma_var);
    part_ginvest_var_isp = var(ginvest);
    part_q_var_isp = var(q);
    part_r_e_var_isp = var(r_e);


    % get simulated data with only the tau shock turned on

    irf_shock = char('eps_tau');

    [zdatapiecewise2_, zdatass_, oo00_,  M00_, violvecbool_last_shock] = ...
        solve_two_constraints_add_violvecbool_no_parse(M00_,M10_,M01_,M11_,oo00_,...
        options00_,constraint1, constraint2,...
        constraint_relax1, constraint_relax2,...
        shockssequence(:,3),irf_shock,nperiods,0,maxiter);

    ytot_pos = strmatch('ytot',M00_.endo_names,'exact');
    gamma_var_pos = strmatch('gamma_var',M00_.endo_names,'exact');


    f0_ytot_data = zdatapiecewise2_(:,ytot_pos);
    f0_ytot_ss = zdatass_(ytot_pos);
    f0_ltot2ytot_data = zdatapiecewise2_(:,ltot2ytot_pos);
    f0_ginvest_data = zdatapiecewise2_(:,ginvest_pos);
    f0_ginvest_ss = zdatass_(ginvest_pos);
    %f0_gamma_var_data = zdatapiecewise2_(:,gamma_var_pos);
    f0_q_data = zdatapiecewise2_(:,q_pos);
    f0_r_e_data = zdatapiecewise2_(:,r_e_pos);


    gdp = 100*bpass(f0_ytot_data/f0_ytot_ss,6,32);
    ltot2ytot = 100*bpass(f0_ltot2ytot_data,6,32);
    %gamma_var = 100*bpass(f0_gamma_var_data,6,32);
    ginvest = 100*bpass(f0_ginvest_data/f0_ginvest_ss,6,32);
    q = 100*bpass(f0_q_data,6,32);
    r_e =100*bpass(f0_r_e_data,6,32);


    part_gdp_var_tau = var(gdp);
    part_ltot2ytot_var_tau = var(ltot2ytot);
    %part_gamma_var_tau = var(gamma_var);
    part_ginvest_var_tau = var(ginvest);
    part_q_var_tau = var(q);
    part_r_e_var_tau = var(r_e);

    %%

    sum_gdp_var = part_gdp_var_a + part_gdp_var_isp + part_gdp_var_tau;
    sum_ltot2ytot_var = part_ltot2ytot_var_a + part_ltot2ytot_var_isp + part_ltot2ytot_var_tau;
    %sum_gamma_var = part_gamma_var_a + part_gamma_var_isp + part_gamma_var_tau;
    sum_ginvest_var = part_ginvest_var_a + part_ginvest_var_isp + part_ginvest_var_tau;
    sum_q_var = part_q_var_a + part_q_var_isp + part_q_var_tau;
    sum_r_e_var = part_r_e_var_a + part_r_e_var_isp + part_r_e_var_tau;


    %%

    gdp_var_share = [round(100*part_gdp_var_a/sum_gdp_var)
        round(100*part_gdp_var_isp/sum_gdp_var)
        100-round(100*part_gdp_var_a/sum_gdp_var)-round(100*part_gdp_var_isp/sum_gdp_var)
        ];

    ltot2ytot_var_share = [round(100*part_ltot2ytot_var_a/sum_ltot2ytot_var)
        round(100*part_ltot2ytot_var_isp/sum_ltot2ytot_var)
        100-round(100*part_ltot2ytot_var_a/sum_ltot2ytot_var)-round(100*part_ltot2ytot_var_isp/sum_ltot2ytot_var)
        ];

    % gamma_var_share =[round(100*part_gamma_var_a/sum_gamma_var)
    %     round(100*part_gamma_var_isp/sum_gamma_var)
    %     100-round(100*part_gamma_var_a/sum_gamma_var)-round(100*part_gamma_var_isp/sum_gamma_var)
    %     ];

    ginvest_var_share =[round(100*part_ginvest_var_a/sum_ginvest_var)
        round(100*part_ginvest_var_isp/sum_ginvest_var)
        100-round(100*part_ginvest_var_a/sum_ginvest_var)-round(100*part_ginvest_var_isp/sum_ginvest_var)
        ];

    q_var_share =[round(100*part_q_var_a/sum_q_var)
        round(100*part_q_var_isp/sum_q_var)
        100-round(100*part_q_var_a/sum_q_var)-round(100*part_q_var_isp/sum_q_var)
        ];


    r_e_var_share =[round(100*part_r_e_var_a/sum_r_e_var)
        round(100*part_r_e_var_isp/sum_r_e_var)
        100-round(100*part_r_e_var_a/sum_r_e_var)-round(100*part_r_e_var_isp/sum_r_e_var)
        ];

    %%

    % texfilename = 'table_calibration';
    % columnlabels = char('AR(1) param.','Innov. St. Dev.');
    % rowlabels = char('TFP','ISP','Volatility$^*$');
    % table_mat = [thetaunc([4,5,6]), thetaunc([1,2,3])];
    % column_precision = [2, 4];
    % maketable(table_mat,texfilename,columnlabels,rowlabels,sideways_switch,caption,column_precision)


    texfilename = 'table_calibration_share';
    columnlabels = char('var(GDP)','var(invest.)','var(invest. p.)','var(credit/GDP)');
    rowlabels = char('TFP','ISP','Volatility');
    table_mat_share =[gdp_var_share, ginvest_var_share, q_var_share, ltot2ytot_var_share];
    column_precision = [0 0 0 0 0];
    maketable(table_mat_share,texfilename,columnlabels,rowlabels,sideways_switch,caption,column_precision)

else
    display(['The model moments cannot be computed'])
end