clear; clc; close all; 
addpath(genpath(pwd));

s = RandStream('mt19937ar','Seed',0);
RandStream.setGlobalStream(s);

[xx,tt] = xlsread(['data_oil.xlsx'],5);
oil_real  = xx(:,1);
oil_supp  = xx(:,2);

%% set the parameters of the exercise
H_grid        = .5;
kernshape     = 'gauss';
% phi_grid      = [10^(-10) 10^(-5) 10^(-4) 10^(-3) 10^(-2):.03:1].^.5;
% optimize_step = 1; %determines how frequently you want to optimize tightness
% H_gridsize    = length(H_grid);
% phi_gridsize  = length(phi_grid);
mods          = 1;
% NN_VAR        = 3;% # di variabili di cui tenere conto nella previsione/scelta del fit

%% prepare the data for the analysis
load data_20_78
if mods==1
        gamma_mean   = prior_MED;
        dataselect   = data_MED(1:T,:);
elseif mods==2
        gamma_mean   = prior_LAR;
        dataselect   = data_LAR(1:T,:);
end
        
% now append the industrial poductiond data
IPI   = {'IPBUSEQ','IPCONGD','IPDCONGD','IPDMAT','IPFINAL','IPMAT','IPNCONGD','IPNMAT'};
disp(IPI)

for jj=1:size(IPI,2);
    pos_IPI(jj) = find(strcmp(IPI{jj},codes_ALL));
end
prior_IPI  =  prior_LAR(pos_IPI,1);
data_IPI   =  data_LAR(1:T,pos_IPI);
gamma_mean = [gamma_mean;prior_IPI];
dataselect = [dataselect data_IPI];
mycodes      = [codes_MED IPI];
m_temp = size(dataselect,2);

% REPLACE NOMINAL WITH REAL OIL
% select 133 onwards (1970 onwards); 
YY_temp                = dataselect(133:end,[14,setdiff(1:m_temp,14)]);
T_temp                 = size(YY_temp,1);
% replace first variable with oil_real
YY_temp(:,1)           = oil_real(1:T_temp,1);
% and add oil production last
YY_lev                 = [YY_temp oil_supp(1:T_temp,1)];
m                      = size(YY_lev,2); % number of variables in the VAR
gamma_mean_d           = zeros(m,1);
codes_ALL_r            = [mycodes(14) mycodes(1,setdiff(1:m_temp,14)) 'OIL_SUPP'];  
lags         = 13;
lev          = [1 3 8 13 14 20];

%% TREAT DATA
YY_tt_d             = YY_lev(2:end,:);
todiff              = setdiff(1:m,[1 3 8 13 14 20]);
for hh=todiff
    YY_tt_d(:,hh)   =  YY_lev(2:end,hh)-YY_lev(1:end-1,hh);
end

d1                     = [diag(gamma_mean_d);zeros(m*(lags-1)+1,m)];
T_est                  = size(YY_tt_d,1)-lags;
[y_l,x_l]                       = matrix_var(YY_tt_d,lags,T_est+lags,m);
[sigu_ar,Omega,OmegaprimeOmega] = get_Omega(T_est,m,lags,y_l,x_l);
phi  = 0.1;
errors = NaN(T_est,size(y_l,2));

for tt=lags+1:T_est+lags
    disp(tt)
    count_tt=tt-lags;
    % WEIGHTS
    H       = T_est.^H_grid;
    Wt      = w_gauss(tt,T_est,H,2);
    big_X_W = x_l'*Wt*x_l;
    big_Y_W = x_l'*Wt*y_l;
    b                      = (big_X_W+OmegaprimeOmega/(phi^2))\(big_Y_W+(OmegaprimeOmega/(phi^2))*d1);
 
    StableDraw = check_stability(b,m,lags);
    if StableDraw ==0
            disp(['at tt=',num2str(tt),' we have an unstable estimate'])
    end
    errors(count_tt,:)     = y_l(count_tt,:)-x_l(count_tt,:)*b;
    slope_b(:,:,count_tt)  = b;
end

% IDENTIFY WITH sign_restrictions
sq               = @(x) squeeze(x);
[Sign_spec,temp] = xlsread('signs_r');
Sign.SIGNS       = Sign_spec(:,1);
Sign.per_SIGNS   = Sign_spec(:,2);

% Choose percentiles for IR and var decomp (keep median first for charts):
vecPerc    = [16 50 85];
vecTime    = 1970.25 : 0.25 : 2014.00;

% details for IRFs with Sign restrictions
Npar = size(slope_b,1);
IRF_hor = 36;
H_sv    = .6;

% estimate time varying volatilities and identify shocks
% Inputs for Sign Restrictions
Sign.h_IRF              = 36;
Sign.max_check_periods  = max(max(Sign.per_SIGNS));
Sign.MaxNumberTry       = 250;
Nshocks=size(Sign.SIGNS,2);
%
KEEP   = 250;
intN   = m;
VCV_TV = [];
count_t=0;
ImpResp_oil_tt = [];
OilDemElast    = [];

% load BVAR_oil_res_delta


%%
disp(strcat('Wait while I am estimating reduced form parameters'));
for tt         = lags+1:24:T_est+lags
    count_t=count_t+1;
    H          = T_est.^H_sv;
    Wt         = w_gauss(tt,T_est,H,2);
    VCV        = errors'*Wt*errors;
    VCV_TV(:,:,tt)= VCV;
    vecBetas   = slope_b(1:end-1,:,tt-lags);
    [ImpResp,A0_id]   = IRFs_oil_largeVAR_signs(vecBetas,VCV,KEEP,intN,lags,Sign);    
    %standardize to have a 1% real increase in the price of oil
    tosave = strcat('ImpResp_oil_',num2str(count_t));
    ImpResp_oil_tt = sq(ImpResp(:,1,:,:));
    eval(['save ',tosave,' ImpResp_oil_tt']);
    clear ImpResp_oil_tt
end

%% PLOT IRFs
ALLIPI   = {'Total','Business Equipment','Consumer Goods','Consumer Durable','Consumer Non-Durable','Final Products','Materials','Durable Materials','Non-Durable Materials'};

clear ImpResp*
close all
ImpResp_oil=[];
m = 29;
Sign.h_IRF = 36;
ImpResp_oil_perc = [];
ImpResp_oil_st = [];
counttt=0;
for tt         = lags+1:24:T_est+lags
    counttt=counttt+1;
    disp(counttt)
    saved = strcat('ImpResp_oil_',num2str(counttt));
    eval(['load ',saved]);
    ImpResp_OIL     = median(squeeze(ImpResp_oil_tt(1,1,:)));
    ImpResp_oil_tt_Stand = ImpResp_oil_tt./ImpResp_OIL; 
    ImpResp_oil_perc_TV(:,:,counttt,:) = cumsum(squeeze(prctile(ImpResp_oil_tt_Stand,[5 50 95], 3)),2);
    ElastDem_TV(counttt,:) = (squeeze(ImpResp_oil_tt_Stand(end,1,:)./ImpResp_oil_tt_Stand(1,1,:)));
end
%
select = [9 21:28];
close all
count=0;
Years =1972:2:2014;
figure
for ll=select
    count=count+1;
    subplot(3,3,count);
    hh=plot(squeeze((ImpResp_oil_perc_TV(ll,1,:,:))),'color','b','linewidth',3);
    hold on
    plot(zeros(1,22),'color','g')
    set(hh([1 3]),'linestyle','--')
    set(gca,'xlim',[1 22])
    set(gca,'ylim',[-.4 0.1])
    set(gca,'XTick',[1:4:22])    
    aa=get(gca);
    title(ALLIPI(count))
    set(gca,'xticklabels',Years(aa.XTick),'Fontsize',18)
end


count=0;
figure
for ll=select
    count=count+1;
    subplot(3,3,count);
    hh=plot(squeeze((ImpResp_oil_perc_TV(ll,12,:,:))),'color','r','linewidth',3)
    hold on
    plot(zeros(1,22),'color','g','linewidth',3)
    set(hh([1 3]),'linestyle','--')
    set(gca,'xlim',[1 22])
    set(gca,'ylim',[-.4 0.1])
    set(gca,'XTick',[1:4:22])    
    aa=get(gca);
    title(ALLIPI(count))
    set(gca,'xticklabels',Years(aa.XTick),'Fontsize',18)
end


%
figure
hh=plot(prctile(ElastDem_TV,[16 50 84],2),'color','r','linewidth',3);
hold on
plot(zeros(1,22),'color','g','linewidth',3)
set(hh([1 3]),'linestyle','--')
set(gca,'xlim',[1 22])
set(gca,'ylim',[-.4 0.1])
set(gca,'XTick',[1:4:22])    
set(gca,'xticklabels',Years(aa.XTick),'Fontsize',18)

