% *****************************************************************************
%  Universite de Lausanne - HEC - MScE - Master thesis - Spring 2014
%  ----------------------------------------------------------------------------
%  DSGE models for developing countries: an application to the Moroccan economy
%  Mohammed Ait Lahcen
% *****************************************************************************

clc; clear all; close all;

setenv('LD_LIBRARY_PATH', '');
addpath('/usr/share/dynare/matlab');


% ----------------------------------------------------------
% Importing data
% ----------------------------------------------------------

% -------- Sample specification -----------------
dates       = [2009.00 : 0.25 : 2019.75];
startdate   = 2009.00; 
enddate     = 2019.75;

startindex  = find(dates == startdate); % position of starting year of sample
endindex    = find(dates == enddate); % position of ending year of sample
xxaxis  = [startdate : 0.25 : enddate]';

% -------- Load the data -------------------------
[data]  = xlsread('LALA',1,'B2:Q45');
T       = size(data,1);
N       = size(data,2);

Y      = data(:,1);    % Nominal GDP
C       = data(:,2);    % Real consumption
I       = data(:,4);    % Real investment
G       = data(:,3);    % Real government expenditures
X       = data(:,5);   % Exports in volume
M       = data(:,6);   % Imports in volume
CPI     = data(:,7);   % CPI inflation
R       = data(:,8);   % Interbank interest rate
Q       = data(:,9);   % Effective exchange rate
E       = data(:,10);   % Nominal exchange rate
C_star  = data(:,12);   % Euro area real GDP
CPI_star  = data(:,14);   % Euro area real GDP
R_star  = data(:,15);   % Euro area short-term interest rates
S  = data(:,16);   % Euro area short-term interest rates

plot(dates,S,'g');
hold off;

% -------- Computing average annual growth ------------------

[growth_data] = [Y C I G X M];
growth = zeros(T,size(growth_data,2));

for j = 1 : size(growth_data,2)
    for jjj = 5 : T
        growth(jjj,j) = log(growth_data(jjj,j)) - log(growth_data(jjj-4,j));
    end
end

mean(growth)*100;
std(growth);
plot(dates,growth(:,1:4));

% -------- Computing long-term shares ------------------

C_Y = C./Y;        % Share of consumption in GDP
mean(C_Y);
plot(dates,C_Y);

G_Y = G./Y;        % Share of government spending in GDP
mean(G_Y);
plot(dates,G_Y);

X_Y = X./Y;        % Share of exports in GDP
mean(X_Y);
M_Y = M./Y;        % Share of imports in GDP
mean(M_Y);

plot(dates,X_Y,'g');
hold all;
plot(dates,M_Y,'r');
hold off;

% -------- Transforming the data ------------------

y = log(Y);
c = log(C);
i = log(I);
g = log(G);
x = log(X);
m = log(M);
pih = log((1+(CPI/100)).^(1/4)); % Transfroming annual inflation into quarterly inflation
r = log((1+(R/100)).^(1/4)); % Transfroming annual interest rates into quarterly interest rates
q = log(Q);
e = log(E);
s = log(S);
c_star = log(C_star);
pi_star = log((1+(CPI_star/100)).^(1/4));
r_star = log((1+(R_star/100)).^(1/4));
s = log(s);

plot(dates,pih,'r');
hold all;
hold off;

%-----------------------------
% HP filtering variables
%-----------------------------

[y_bar,y_hat] = hpfilter(y,1600);
[c_bar,c_hat] = hpfilter(c,1600);
[i_bar,i_hat] = hpfilter(i,1600);
[g_bar,g_hat] = hpfilter(g,1600);
[x_bar,x_hat] = hpfilter(x,1600);
[m_bar,m_hat] = hpfilter(m,1600);
[pih_bar,pih_hat] = hpfilter(pih,1600);
[r_bar,r_hat] = hpfilter(r,1600);
[q_bar,q_hat] = hpfilter(q,1600);
[e_bar,e_hat] = hpfilter(e,1600);
[s_bar,s_hat] = hpfilter(s,1600);
[c_star_bar,c_star_hat] = hpfilter(c_star,1600);
[pi_star_bar,pi_star_hat] = hpfilter(pi_star,1600);
[r_star_bar,r_star_hat] = hpfilter(r_star,1600);

plot(dates,pih_hat,'r');


%-----------------------------
% Business cycle figures
%-----------------------------

vlist   = {'y','c','x','m'};
vnames  = {'Total Output','Consumption','Exports','Imports'};

figure(1)
for j = 1 : length(vnames)
    subplot(2,2,j)
    [AX,H1,H2] = plotyy(dates,eval(vlist{j}),dates,eval(strcat(vlist{j},'_hat')));
    H3 = line(dates,eval(strcat(vlist{j},'_bar')),'Parent',AX(1));
    H4 = line(dates,0,'Parent',AX(2));
    set(AX,'xlim',[startdate enddate],'xtick',[startdate:2:enddate]);
    set(AX(1),'ylim',[min(eval(vlist{j}))-0.1 max(eval(vlist{j}))+0.1]);
    set(AX(2),'ylim',[-0.3,0.3]);
    set(AX(1),'YTick',[1:0.25:15]);
    set(AX(2),'YTick',[-0.3:0.1:0.3])
    set(H1,'Color','b');
    set(H2,'LineStyle','-','Color','r');
    set(H3,'Color','g');
    title(vnames{j});
    axis manual
    %legend('Log','HP filtered trend','HP filtered cycle');
    hold all
end


vlist   = {'i','g','x','m'};
vnames  = {'Investment','Government spending','Exports','Imports'};

figure(2)
for j = 1 : length(vnames)
    subplot(2,2,j)
    [AX,H1,H2] = plotyy(dates,eval(vlist{j}),dates,eval(strcat(vlist{j},'_hat')));
    H3 = line(dates,eval(strcat(vlist{j},'_bar')),'Parent',AX(1));
    H4 = line(dates,0,'Parent',AX(2));
    set(AX,'xlim',[startdate enddate],'xtick',[startdate:2:enddate]);
    set(AX(1),'ylim',[min(eval(vlist{j}))-0.1 max(eval(vlist{j}))+0.1]);
    set(AX(2),'ylim',[-0.3,0.3]);
    set(AX(1),'YTick',[1:0.25:15]);
    set(AX(2),'YTick',[-0.3:0.1:0.3])
    set(H1,'Color','b');
    set(H2,'LineStyle','-','Color','r');
    set(H3,'Color','g');
    title(vnames{j});
    axis manual
    %legend('Log','HP filtered trend','HP filtered cycle');
    hold all
end


vlist   = {'pih','pi_star','r','e','q','s'};
vnames  = {'CPI inflation','CPI foreign inflation','Nominal interest rate','Nominal exchange rate','Real exchange rate','Terms of trade'};

figure(3)
for j = 1 : length(vnames)
    subplot(3,2,j)
    AX = plot(dates,eval(vlist{j}),'b',dates,eval(strcat(vlist{j},'_bar')),'-r');
    set(gca,'xlim',[startdate enddate],'xtick',[startdate:2:enddate]);
    %set(gca,'ylim',[min(eval(vlist{j})) max(eval(vlist{j}))]);
    %set(gca,'YTick',[-5:0.1:5])
    title(vnames{j});
    %axis manual
    %legend('Log','HP filtered trend','FontSize',4);
    hold all
end
hold off



%-----------------------------
% Business cycle statistics
%-----------------------------

[data_stat] = [y_hat c_hat i_hat g_hat x_hat m_hat pih_hat r_hat q_hat e_hat s_hat];

V               = size(data_stat,2);
[data_acf]      = zeros(5,V);
[data_adf]      = zeros(1,V);
[data_std]      = zeros(1,V);
[data_relstd]   = zeros(1,V);
[data_corr]     = zeros(1,V);

for t = 1:V
data_acf(:,t)   = autocorr(data_stat(:,t),4);
data_adf(1,t)   = adftest(data_stat(:,t),'model','ARD');  % ADF test for stationarity
data_std(1,t)   = std(data_stat(:,t));
data_relstd(1,t)= std(data_stat(:,t))/std(data_stat(:,1));
data_corr(1,t)  = corr(data_stat(:,1),data_stat(:,t));
end

data_std    = round(data_std*10000)/100 ;
data_relstd = round(data_relstd*100)/100 ;
data_acf    = round(data_acf*100)/100 ;
data_corr   = round(data_corr*100)/100 ;

disp('Y C I G X M R Q E S')
disp(data_acf)


%-----------------------------------
% Euro area AR(1) processes
%-----------------------------------

%model = arima('ARLags',1,'Constant',0);
%c_star_ar1 = estimate(model,c_star_hat)
%pi_star_ar1 = estimate(model,pi_star_hat)
%r_star_ar1 = estimate(model,r_star_hat)


%-----------------------------------
% Storing observables for estimation
%-----------------------------------

[data_est] = [y c pih s q];
save('data_est.mat', 'data_est')

dynare bol.mod
%%%%%%%%%%%%%
%-----------------------------------
% Computing the theoretical moments
%-----------------------------------

load('bol_results.mat');

vlist   = {'y','c','pih','r','q','s'};
for j = 1 : length(vlist)
model_acf(:,j)   = autocorr(eval(strcat('oo_.SmoothedVariables.Mean.',vlist{j})),4);
model_adf(1,j)   = adftest(eval(strcat('oo_.SmoothedVariables.Mean.',vlist{j})),'model','ARD');  % ADF test for stationarity
model_std(1,j)   = std(eval(strcat('oo_.SmoothedVariables.Mean.',vlist{j})));
model_relstd(1,j)= std(eval(strcat('oo_.SmoothedVariables.Mean.',vlist{j})))/std(eval(strcat('oo_.SmoothedVariables.Mean.',vlist{1})));
model_corr(1,j)  = corr(eval(strcat('oo_.SmoothedVariables.Mean.',vlist{1})),eval(strcat('oo_.SmoothedVariables.Mean.',vlist{j})));
end


%-----------------------------------
% Plotting historical simulation
%-----------------------------------

load('bol_results.mat');


vlist   = {'y','c','pih','r','q','s'};
vnames  = {'Total Output','Consumption','Inflation','Interest rate','Real exchange rate','Terms of trade'};

% % Using updated variables
% figure('Name','Using updated variables')
% for j = 1 : length(vnames)
%     subplot(3,3,j)
%     AX = plot(dates,eval(strcat(vlist{j},'_hat')),'b',dates,eval(strcat('[0; oo_.UpdatedVariables.Mean.',vlist{j},']')),'-r');
%     set(gca,'xlim',[startdate enddate],'xtick',[startdate:2:enddate]);
%     title(vnames{j});
%     hold all
% end
% hold off
% 
% % Using filtered variables (one step ahead forecast)
% figure('Name','Using filtered variables')
% for j = 1 : length(vnames)
%     subplot(3,3,j)
%     AX = plot(dates,eval(strcat(vlist{j},'_hat')),'b',dates,eval(strcat('[0; oo_.FilteredVariables.Mean.',vlist{j},']')),'-r');
%     set(gca,'xlim',[startdate enddate],'xtick',[startdate:2:enddate]);
%     title(vnames{j});
%     hold all
% end
% hold off

% Using smoothed variables (using all information)
figure('Name','Using smoothed variables (1)')
for j = 1 : 4
    subplot(2,2,j)
    AX = plot(dates,eval(strcat(vlist{j},'_hat')),'b',dates,eval(strcat('[0; oo_.SmoothedVariables.Mean.',vlist{j},']')),'--r');
    set(gca,'xlim',[startdate enddate],'xtick',[startdate:2:enddate]);
    title(vnames{j},'fontsize',16);
    hold all
end
hold off

%-----------------------------------
% Plotting Bayesian IRFs
%-----------------------------------

% dynare NK_SOE_model_v9_sim.mod noclearall                         % Launch simulation using new posterior parameters
load('bol_results.mat');                % Load simulation results

sim_irf = oo_.irfs;

load('bol_results.mat');                    % Loading Bayesian estimation results

post_param = oo_.posterior_mean.parameters;             % Loading posterior parameters
save('post_param.mat', '-struct','post_param');                   % Save posterior parameters for simulation

post_shock = oo_.posterior_mean.shocks_std;             % Loading posterior shock std

[post_shock.sigma_at] = post_shock.eps_at;              % Renaming shocks to load in Dynare
post_shock = rmfield(post_shock,'eps_at');
[post_shock.sigma_an] = post_shock.eps_an;
post_shock = rmfield(post_shock,'eps_an');
[post_shock.sigma_r] = post_shock.eps_r;
post_shock = rmfield(post_shock,'eps_r');
save('post_shock.mat', '-struct','post_shock');                   % Save posterior shock std for simulation

post_irf = oo_.PosteriorIRF.dsge;

hor = 40;
xx = zeros(hor,1);


slist   = {'eps_at','eps_an','eps_r','eps_cstar','eps_pistar'};
snames  = {'Tradables productivity shock','Non-tradables productivity shock','Interest rate shock','World consumption shock','Imported inflation shock'};

vlist   = {'y','y_t','y_nf','y_ni','rw_f','rw_i','c','c_th','c_tf','pi','pi_t','pi_nf','pi_ni','r','q','s'};
vnames  = {'Total Output','Tradables Output','Formal Non-tradables Output/Consumption','Informal Non-tradables Output/Consumption','Real formal wage','Real informal wage','Consumption','Domestic tradables consumption','Imports consumption','CPI inflation','Tradables inflation','Formal non-tradables inflation','Informal non-tradables inflation','Nominal interest rate','Real exchange rate','Terms of trade'};


for j = 1 : length(slist)

post_irf.Mean.(strcat('rw_f_',slist{j})) = eval(strcat('(post_irf.Mean.w_f_',slist{j},')')) - eval(strcat('(post_irf.Mean.p_',slist{j},')'));

post_irf.Mean.(strcat('rw_i_',slist{j})) = eval(strcat('(post_irf.Mean.w_i_',slist{j},')')) - eval(strcat('(post_irf.Mean.p_',slist{j},')'));

post_irf.HPDinf.(strcat('rw_f_',slist{j})) = eval(strcat('(post_irf.HPDinf.w_f_',slist{j},')')) - eval(strcat('(post_irf.HPDinf.p_',slist{j},')'));

post_irf.HPDsup.(strcat('rw_f_',slist{j})) = eval(strcat('(post_irf.HPDsup.w_f_',slist{j},')')) - eval(strcat('(post_irf.HPDsup.p_',slist{j},')'));

post_irf.HPDinf.(strcat('rw_i_',slist{j})) = eval(strcat('(post_irf.HPDinf.w_i_',slist{j},')')) - eval(strcat('(post_irf.HPDinf.p_',slist{j},')'));

post_irf.HPDsup.(strcat('rw_i_',slist{j})) = eval(strcat('(post_irf.HPDsup.w_i_',slist{j},')')) - eval(strcat('(post_irf.HPDsup.p_',slist{j},')'));

sim_irf.(strcat('rw_f_',slist{j})) = eval(strcat('(sim_irf.w_f_',slist{j},')')) - eval(strcat('(sim_irf.p_',slist{j},')'));

sim_irf.(strcat('rw_i_',slist{j})) = eval(strcat('(sim_irf.w_i_',slist{j},')')) - eval(strcat('(sim_irf.p_',slist{j},')'));

end



    for j = 1 : length(snames)

         
        
    for jj = 1 : 4
        
     
    figure('name',strcat(snames{j},'-',num2str(jj)))
         
            
    for jjj = (1+(length(vnames)/4*(jj-1))) : ((length(vnames)/4)+(length(vnames)/4*(jj-1)))

        scalea = min(eval(strcat('post_irf.HPDinf.',vlist{jjj},'_',slist{j}))).*100;

        subplot(2,2,jjj-(length(vnames)/4*(jj-1)))
            area(0:1:hor-1,eval(strcat('post_irf.HPDsup.',vlist{jjj},'_',slist{j})).*100,scalea,'FaceColor',[.7 .7 .7],...
                'EdgeColor','k','LineWidth',.5,'Linestyle','none');
        hold on
            area(0:1:hor-1,eval(strcat('post_irf.HPDinf.',vlist{jjj},'_',slist{j})).*100,scalea,'FaceColor',[1 1 1],...
                'EdgeColor','k','LineWidth',.5,'Linestyle','none');
        hold on

        % -----------------------------------------------
        % The main figure
        % -----------------------------------------------
        jm = plot(0:1:hor-1,eval(strcat('post_irf.Mean.',vlist{jjj},'_',slist{j})).*100,'r');
        %set(jm(1),'linewidth',1,'Color','r')
       
        jmm = plot(0:1:hor-1,eval(strcat('sim_irf.',vlist{jjj},'_',slist{j})).*100,'b');
        %set(jm(1),'linewidth',1)

        jxx = plot(0:1:hor-1,xx','k');
        set(jxx(1),'linewidth',0.5)
                
        title(vnames{jjj},'fontsize',12)
        xlabel('Quarters after shock','fontsize',10)
        ylabel('Percent deviation','fontsize',10)
        
        %legend show;

        axis tight
        set(gca,'Layer','top')
        
        
        
    end
    
    print('-dpng',strcat(snames{j},'-',num2str(jj)))
    
    end

    end

    % For legend see http://stackoverflow.com/questions/13685967/how-to-show-legend-for-only-a-specific-subset-of-curves-in-the-plotting
    
    
