%% SOE model with DNWR
close all
clearvars
clc

% restoredefaultpath
% warning('off')
% addpath(genpath('routines'))
% addpath(genpath('C:/Dynare/5.4'))
% addpath(genpath('/Applications/Dynare/5.4'))

setup = 'IRF_SOE';

inLevels    = 1;   % whether to run in levels or in logs
useOccBin   = 1;

if inLevels==0
    eval(['dynare SOE'])
elseif inLevels==1
    eval(['dynare SOE_levels'])
end

N   = length(M_.endo_names);     % number of endogenous variables
pos = M_.endo_names;   % position of endogenous variables

% Structure with parameter values
par.alpha   = 0.35;         % capital share
par.beta    = 1/(1.005);       % discount factor 
par.sig     = 1.5;           % IES
par.eta     = .999;        % DNWR

par.iota    = 1e-4;         % penalty cost for bond holdings

% par.rho     = 0.99;         % shock persistence
% par.rho2    = 0;
par.rho     = 1.64;
par.rho2    = -0.65; 

par.om      = 0.7;         % home bias
par.psi     = 2.5;         % Armington elasticity

par.noMigr  = 1;           % 1: no migration case
par.gam     = 0.2;         % migration elasticity       
par.ss_n    = 0.01;        % share of migrants (adjusted for share_w below)

par.occbin_DNWR_bind = 0;

% Implied parameters
par.share_w = 1-par.alpha; % share of workers in total population
par.ss_n    = par.ss_n/par.share_w; 
par.ss_Nk   = par.alpha;
par.ss_K    = 1;

% feed parameters into dynare which saves parameters in M_.params
for i = 1:length(M_.param_names)
    M_.params(i) = par.(M_.param_names{i});  % some parameter values saved in 'par'
end

% variance of shock
M_.Sigma_e(1,1) = (0.005)^2; % variance of shock (IRF: one-std-dev shock)

% if estimating in levels, feed in steady state
if inLevels==1
    % Structure with steady-state values
    ss.pQ   = 1;
    ss.e    = 1;
    ss.Q    = 1;
    ss.L    = 1;
    ss.N    = 1;
    ss.Nw   = 1-par.alpha;
    ss.Z    = par.alpha^(-par.alpha) * (1-par.alpha)^(par.alpha-1);
    ss.r    = 1;
    ss.w    = 1;
    ss.C    = ss.Q;
    ss.cw   = 1;
    ss.ck   = 1;
    ss.M    = 1-par.om;
    ss.X    = ss.M;
    ss.D    = par.om;
    ss.B    = 0;
    ss.V    = 0;
    ss.w_lag= ss.w;
    ss.e_lag= ss.e;
    ss.Z_lag= ss.Z;

    for i = 1:N
        oo_.steady_state(i) = ss.(M_.endo_names{i});
    end
    % adjust variance of shock s.t. 1 std dev = 0.01
    M_.Sigma_e = ss.Z^2 * M_.Sigma_e;
end

if useOccBin==1
    nperiods = 100;
    randn('seed',3);
    shockssequence = 1*randn(nperiods,1)*sqrt(M_.Sigma_e);
    for i = 1:nperiods
        M_.surprise_shocks(i).periods = i;
        M_.surprise_shocks(i).value = shockssequence(i);
        M_.surprise_shocks(i).exo_id  = 1;
    end

    % get policy function: oo_.dr
    options_occbin_.maxit = 200;
    options_occbin_.simul_check_ahead_periods=500;
    [M_,options_] = occbin.setup(M_,options_, options_occbin_);
%     options_.occbin.simul.periodic_solution = 1;
    oo_ = occbin.solver(M_, oo_, options_);
    for i = 1:N
        irf.(pos{i})(1,:) = 100*(oo_.occbin.simul.piecewise(:,i)-oo_.steady_state(i));
    end
else
    % get policy function: oo_.dr
    options_.order  = 1; options_.nograph= 1; options_.irf = 40;
    [~, oo_, options_] = stoch_simul(M_, options_, oo_,M_.endo_names);
    for i = 1:N
        irf.(pos{i})(1,:) = 100*oo_.irfs.([pos{i},'_eps_Z']);
    end
end

% save irf: oo_.irf
for i = 1:N
    % if in logs, irfs are in relative deviations
    % if in levels, irfs are in absolute deviations -> transform into
    % relative deviations
    if inLevels==1 & oo_.steady_state(i)~=0
        irf.(pos{i})(1,:) = irf.(pos{i})(1,:) / oo_.steady_state(i);
    end
    irf.(pos{i})(1,abs(irf.(pos{i})(1,:))<1e-10) = 0;  % set small values to 0
end

irf.LN = irf.L + irf.Nw;
irf.W = irf.w + irf.e;


%% Plot IRF
hFig=figure(1);
set(gcf, 'Color', 'w'); % set background color to white
set(hFig,'Units', 'normalized','Position', [0 0 .8 .8]) % size of figure

% Choose which variables to plot + add variable names
var = {'Z','C','L','e','pQ','w','N','B','W'};
s1 = 3; s2 = 3;

% var = {'Z','L','N','W'};
% s1 = 2; s2 = 2;


varName.C   = 'Consumption';
varName.L   = 'Labor per capita';
varName.e   = 'Real exchange rate';
varName.pQ  = 'ToT';
varName.w   = 'Real wage';
varName.W   = 'Nominal wage';
varName.N   = 'Population';
varName.B   = 'Bond holdings';
varName.Z   = 'TFP';
varName.LN  = 'Total labor';


for j = 1:length(var)
    v = var{j};     % variable to be plotted
    T = length(irf.(v));    % length of IRFs

    subplot(s1,s2,j)  % 2x4 subplots
    hold on     % to plot several lines 
    plot(1:T,irf.(v),'Linewidth',3.5)
    plot(1:T,zeros(1,T),'k--','Linewidth',2)    % add reference line
    hold off
    set(gca, 'Fontsize', 16); % default font size
    grid minor  % add minor grid lines 
    axis tight  % avoid gaps between line and axes
    xlabel('quarter','FontSize',20)
    if strcmp(v,'infl') | strcmp(v,'int') | strcmp(v,'r')
        ylabel('percentage points','FontSize',20)
    else
        ylabel('percent','FontSize',20)
    end
    title(varName.(v),'FontSize',26)    
end

[abs(mean(irf.L)) std(irf.L) abs(mean(irf.L))/std(irf.L)]

