% Gali-van Rens (2010) model
% TvR, 12/2008, revised 05/2010


clear;
close all;

addpath ..\..\dynare\4.0.3\matlab    % for code to run on server


%----------------------------------------------------------------
% Settings
%----------------------------------------------------------------

HP=0; % Switch on for HP filtering (lambda=1600)
% First (or second) order approximation, set in mod-file
% Number of periods for simulations, set in mod-file
% 1000 initial periods dropped, set below


%----------------------------------------------------------------
% Solve the model
%----------------------------------------------------------------

Iter = 0;
StartTime = cputime;
results = [];
save resultstable results;

% Calibration:

make_plots = 0;

Rbar_vals = [0 .95];
for Rbar_index = 1:size(Rbar_vals,2)
Rbar = Rbar_vals(Rbar_index);
Rbar_bak = Rbar;
Rcurv = 1;     % needs to be non-negative integer: 0 is flex, 1 is quadratic

%psi_vals = .1;
psi_vals = 0.3; %[.05 .1 .2 .3];
for psi_index = 1:size(psi_vals,2)
psi = psi_vals(psi_index);

Nstst_vals = 0.7;
for Nstst_index = 1:size(Nstst_vals,2)
Nstst = Nstst_vals(Nstst_index);

mu_vals = 1; %[1 2]; %0.6; 
%mu_vals = [0.2 0.6 0.8];
for mu_index = 1:size(mu_vals,2)
mu = mu_vals(mu_index);

if Rbar==0.95
    sigmaa_vals = 0.00186;
elseif Rbar==0
    sigmaa_vals = 0.00202;
end
for sigmaa_index = 1:size(sigmaa_vals,2)
sigmaa = sigmaa_vals(sigmaa_index);

%sigmaz_vals = 1.1*.001828*[.9 1 1.1];
if Rbar==0.95
    sigmaz_vals = .00173;
elseif Rbar==0
    sigmaz_vals = .00208;
end
for sigmaz_index = 1:size(sigmaz_vals,2)
sigmaz = sigmaz_vals(sigmaz_index);

kappa_vals = [1.35 0.816 0.373 0]; %[1.405 0.841 0.380 0.182 0]; % [3.01 2.53 2.11 1.74 1.405 1.108 0.841 0.598 0.380 0.182];
if mu==2
    kappa_vals = 11*kappa_vals;
end
for kappa_index = 1:size(kappa_vals,2)
kappa = kappa_vals(kappa_index);
if kappa == 0
    Rbar = 0;
else
    Rbar = Rbar_bak;
end

save parameterfile psi mu sigmaa sigmaz Rbar Rcurv Nstst kappa;

% Call Dynare to solve the model:
tic;
dynare model_nocapital noclearall;
toc;


%----------------------------------------------------------------
% Table business cycle stats
%----------------------------------------------------------------

data = [p y n w w_LB w_UB R e v c sdf a z];
data = data(1001:end-3,:);
varnames = ['p   ';'y   ';'n   ';'w   ';'w_LB';'w_UB';'R   ';'e   ';'v   ';'c   ';'sdf ';'a   ';'z   '];
varnames_dyn = M_.endo_names;
if varnames_dyn ~= varnames
    error('you changed the number or order of variables');
end

index_y = find(varnames(:,1)=='y'&varnames(:,2)==' ');
index_n = find(varnames(:,1)=='n'&varnames(:,2)==' ');
index_p = find(varnames(:,1)=='p'&varnames(:,2)==' ');
index_w = find(varnames(:,1)=='w'&varnames(:,2)==' ');
index_wUB = find(varnames(:,1)=='w'&varnames(:,3)=='U');
index_wLB = find(varnames(:,1)=='w'&varnames(:,3)=='L');
index_e = find(varnames(:,1)=='e'&varnames(:,2)==' ');
index_c = find(varnames(:,1)=='c'&varnames(:,2)==' ');
index_v = find(varnames(:,1)=='v'&varnames(:,2)==' ');
index_nolog = find( (varnames(:,1)=='R'&varnames(:,2)==' ') );

levelstst = exp(oo_.steady_state);
levelstst(index_nolog) = log(levelstst(index_nolog));

% THEORETICAL MOMENTS (works only if no variables are dropped by Dynare):
% levelmeans = exp(oo_.mean);
% levelmeans(index_nolog) = log(levelmeans(index_nolog));
% sd = sqrt(diag(oo_.var));
% relsd = sqrt(diag(oo_.var))/sqrt(oo_.var(index_y,index_y));
% autocorr = diag(oo_.autocorr{1,1});
% corry = oo_.var(:,index_y)./(sd*sd(index_y));
% corrn = oo_.var(:,index_n)./(sd*sd(index_n));
% table('THEORETICAL BUSINESS CYCLE MOMENTS)',strvcat('VARIABLE','ST.STATE','MEAN','AUTOCORR','STD.DEV.','REL.S.D.','CORR.Y','CORR.N'),lgy_,[levelstst,levelmeans,autocorr,sd,relsd,corry,corrn],10,8,4);

mean_simul = exp(mean(data))';
mean_simul(index_nolog) = log(mean_simul(index_nolog));

% Filter the data
if HP==1
    for j = 1:size(data,2)
        data(:,j) = data(:,j) - hpfilter(data(:,j),1600); % idealfilter(timeseries(data(:,j)),[1/32 1/6],'pass')
    end
end

sd_simul = sqrt(var(data))';
relsd_simul = sd_simul/sd_simul(index_y);

for j = 1:size(data,2)
    if sd_simul(j)<1e-12
        autocorr_simul(j,1) = NaN;
        corry_simul(j,1) = 0;
        corrn_simul(j,1) = 0;
    else
        autocorr_simul(j,1) = corr(data(2:end,j),data(1:end-1,j));
        corry_simul(j,1) = corr(data(:,j),data(:,index_y));
        corrn_simul(j,1) = corr(data(:,j),data(:,index_n));
    end
end

% Calibration targets:
Eff_stst_fl = ( (psi/(1+phi))*(1/zeta)*(1/(1-PsiF-PsiH)) )^(1/(1+phi));
Nstst_fl = (1-alpha)*(1-PsiF-PsiH)*((1+zeta)/gamma)*levelstst(index_y)^(1-eta);
VacCosts = ((kappa/(1+mu))*(levelstst(index_v)^(1+mu)))/levelstst(index_y);

% Output the results:
disp(sprintf('\nCALIBRATION:'));
disp(sprintf('\nSteady state:'));
disp(sprintf('zeta   = %g ==> Eff_stst_fl = %g (target = 1 in frictionless model)', zeta, Eff_stst_fl));
disp(sprintf('gamma  = %g ==> Nstst_fl = %g (target = %g in frictionless model)', gamma, Nstst_fl, Nstst));
disp(sprintf('Vacancy posting costs as fraction of output = %g', VacCosts));
disp(sprintf('\nSecond moments (for HP filter):'));
disp(sprintf('sigmaa = %g ==> corr(p,y)   = %g pre 1984 (target = 0.77)', sigmaa, corry_simul(index_p)));
disp(sprintf('sigmaz = %g ==> sd(n)/sd(y) = %g pre 1984 (target = 0.66)', sigmaz, relsd_simul(index_n)));
disp(sprintf('psi    = %g ==> sd(n)/sd(y) = %g post 1984 (target = 0.97)', psi, relsd_simul(index_n)));

table('BUSINESS CYCLE MOMENTS (logs ~R, no filter, 50000+1000 simul)',strvcat('VARIABLE','ST.STATE','MEAN','AUTOCORR','STD.DEV.','REL.S.D.','CORR.Y','CORR.N'),varnames,[levelstst,mean_simul,autocorr_simul,sd_simul,relsd_simul,corry_simul,corrn_simul],10,8,4);

% Save the most important results:
load resultstable;
results = [ results ; [psi mu Rcurv Rbar 100*sigmaa 100*sigmaz levelstst(index_n) VacCosts corry_simul(index_p) corrn_simul(index_p) relsd_simul(index_n) relsd_simul(index_w) 100*sd_simul(index_y)] ];
save resultstable results;

% Some plots:

if make_plots

t=1:size(data,1);
period = 1:500;
w_plot = data(period,find(varnames(:,1)=='w'&varnames(:,2)==' '));
w_LB_plot = data(period,find(varnames(:,1)=='w'&varnames(:,3)=='L'));
w_UB_plot = data(period,find(varnames(:,1)=='w'&varnames(:,3)=='U'));

% Accuracy check endogenous wage rigidity
what_plot = (2*w_plot-w_LB_plot-w_UB_plot)./(w_UB_plot-w_LB_plot);
R_nonlin_plot = Rbar*(1-what_plot.^(2*Rcurv));
R_simul_plot = data(period,find(varnames(:,1)=='R'&varnames(:,2)==' '));
figure(1);
plot(what_plot,R_simul_plot,'bo',what_plot,R_nonlin_plot,'rd');
axis([-.8 .8 .8 1.05]);

% Simulated wage and bargaining set
figure(2);
plot(t(period),w_plot,'b-',t(period),w_LB_plot,'r--',t(period),w_UB_plot,'r-');

end  % if make_plots

% End loop over parameter values:
Iter = Iter + 1;
disp(sprintf('Iterations completed: %g, Elapsed time: %g',Iter,cputime-StartTime));
end
end
end
end
end
end
end


%----------------------------------------------------------------
% Output results table
%----------------------------------------------------------------

load resultstable;
modelnr = int2str([1:size(results,1)]');
table('RESULTS FOR VARIOUS CALIBRATIONS',strvcat('MODEL','PSI','MU','RCURV','RBAR','SIGMAA','SIGMAZ','N_STST','VACCOSTS','CORR_PY','CORR_PN','RELSD_N','RELSD_W','SD_Y'),modelnr,results,10,8,4);