clear;clc;close all;

%-------------------------------------------------------------------
% call getdata.m to fetch the data, adjust the getdata.m function if
% different / more variables shall be used
%-------------------------------------------------------------------
[years,y,g,t,r,ex,p,infl,ly,lg,lt,ls,ld,lex,lp,countries,ifs,year,crisis,lrate,mgov,mpers,Dmgov,wardummy,mgov2,exregime,eaunion]=getdata_jst_panel;

%Specify time period
Myear=reshape(year,[years,countries]);
t1=find(Myear(:,1)==1917);              %Start date
t2=find(Myear(:,1)==2016);              %End date

y=cutsample(t1,t2,y);
g=cutsample(t1,t2,g);
t=cutsample(t1,t2,t);
r=cutsample(t1,t2,r);
ex=cutsample(t1,t2,ex);
p=cutsample(t1,t2,p);
infl=cutsample(t1,t2,infl);
ly=cutsample(t1,t2,ly);
lg=cutsample(t1,t2,lg);
mgov2=cutsample(t1,t2,mgov2.*100);
lmg=log(mgov2);
Dmgov=cutsample(t1,t2,Dmgov.*100);
lt=cutsample(t1,t2,lt);
ls=cutsample(t1,t2,ls);
ld=cutsample(t1,t2,ld);
lex=cutsample(t1,t2,lex);
lp=cutsample(t1,t2,lp);
year=cutsample(t1,t2,year);
crisis=cutsample(t1,t2,crisis);
lrate=cutsample(t1,t2,lrate);
ifs=cutsample(t1,t2,ifs);
wardummy=cutsample(t1,t2,wardummy);
exregime=cutsample(t1,t2,exregime);
eaunion=cutsample(t1,t2,eaunion);

years=size(y,1)./countries;

ci = 1.645;  %90%
ci1 = 1;      %68%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Calculate Indicator Variable
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Mr=reshape(r,[years,countries]);
MIr=zeros(size(Mr));

for i=1:countries
   for j=1:size(MIr,1)-1
       %if Mr(j,i)<=meanr(:,i)-1*stdr(:,i)
       %if Mr(j,i)<=quan(:,i)
       if Mr(j,i)<=1
       %if Mr(j,i)<=0
        MIr(j+1,i)=1;
        %elseif Mr(j,i)>meanr(:,i)-1*stdr(:,i)
        %elseif Mr(j,i)>quan(:,i)
        elseif Mr(j,i)>1
        %elseif Mr(j,i)>=0
       MIr(j+1,i)=0;
        end
   end   
end

Mexregime=reshape(exregime,[years,countries]);
MIexregime=zeros(size(Mexregime));
Meaunion=reshape(eaunion,[years,countries]);

for i=1:countries
   for j=1:size(MIexregime,1)-1
            if Meaunion(j,i)==1          %EU Currency Union
        MIexregime(j+1,i)=1;
        end
   end   
end

Ir=reshape(MIr,[countries*years,1]);
Iexregime=reshape(MIexregime,[countries*years,1]);
I=Iexregime;
Iexregime=Iexregime;                                    
                                                        
                                                        
nlags=4; %Number of lags we are interested in +1

MI = reshape(I,[years,countries]);
MI(1:nlags-1,:) =[];
I= reshape(MI,[countries*years-countries*(nlags-1),1]);

MLy=reshape(ly,[years,countries]);
MLy=MLy(nlags-1:end-1,:);
Ly = reshape(MLy,[countries*years-countries*(nlags-1),1]);

ML2y=reshape(ly,[years,countries]);
ML2y=ML2y(nlags-2:end-2,:);
L2y = reshape(ML2y,[countries*years-countries*(nlags-1),1]);

MLg = reshape(lg,[years,countries]);
MLg=MLg(nlags-1:end-1,:);
Lg = reshape(MLg,[countries*years-countries*(nlags-1),1]);

ML2g=reshape(lg,[years,countries]);
ML2g=ML2g(nlags-2:end-2,:);
L2g = reshape(ML2g,[countries*years-countries*(nlags-1),1]);

MLmg = reshape(lmg,[years,countries]);
MLmg=MLmg(nlags-1:end-1,:);
Lmg = reshape(MLmg,[countries*years-countries*(nlags-1),1]);

ML2mg=reshape(lmg,[years,countries]);
ML2mg=ML2mg(nlags-2:end-2,:);
L2mg = reshape(ML2mg,[countries*years-countries*(nlags-1),1]);


MLr = reshape(r,[years,countries]);
MLr=MLr(nlags-1:end-1,:);
Lr = reshape(MLr,[countries*years-countries*(nlags-1),1]);

ML2r=reshape(r,[years,countries]);
ML2r=ML2r(nlags-2:end-2,:);
L2r = reshape(ML2r,[countries*years-countries*(nlags-1),1]);


MLrlong = reshape(lrate,[years,countries]);
MLrlong=MLrlong(nlags-1:end-1,:);
Lrlong = reshape(MLrlong,[countries*years-countries*(nlags-1),1]);

ML2rlong=reshape(lrate,[years,countries]);
ML2rlong=ML2rlong(nlags-2:end-2,:);
L2rlong = reshape(ML2rlong,[countries*years-countries*(nlags-1),1]);

MLex = reshape(lex,[years,countries]);
MLex=MLex(nlags-1:end-1,:);
Lex = reshape(MLex,[countries*years-countries*(nlags-1),1]);

ML2ex=reshape(lex,[years,countries]);
ML2ex=ML2ex(nlags-2:end-2,:);
L2ex = reshape(ML2ex,[countries*years-countries*(nlags-1),1]);

MLp = reshape(lp,[years,countries]);
MLp=MLp(nlags-1:end-1,:);
Lp = reshape(MLp,[countries*years-countries*(nlags-1),1]);

ML2p=reshape(lp,[years,countries]);
ML2p=ML2p(nlags-2:end-2,:);
L2p = reshape(ML2p,[countries*years-countries*(nlags-1),1]);

MLinfl = reshape(infl,[years,countries]);
MLinfl=MLinfl(nlags-1:end-1,:);
Linfl = reshape(MLinfl,[countries*years-countries*(nlags-1),1]);

ML2infl=reshape(infl,[years,countries]);
ML2infl=ML2infl(nlags-2:end-2,:);
L2infl = reshape(ML2infl,[countries*years-countries*(nlags-1),1]);

MLs = reshape(ls,[years,countries]);
MLs=MLs(nlags-1:end-1,:);
Ls = reshape(MLs,[countries*years-countries*(nlags-1),1]);

ML2s=reshape(ls,[years,countries]);
ML2s=ML2s(nlags-2:end-2,:);
L2s = reshape(ML2s,[countries*years-countries*(nlags-1),1]);

Mcountry = reshape(ifs,[years,countries]);
Mcountry(1:nlags-1,:) =[];
country = reshape(Mcountry,[countries*years-countries*(nlags-1),1]);

Mwar = reshape(wardummy,[years,countries]);
Mwar(1:nlags-1,:) =[];
Iwar = reshape(Mwar,[countries*years-countries*(nlags-1),1]);

Mcrisis = reshape(crisis,[years,countries]);
Mcrisis(1:nlags-1,:) =[];
Icrisis = reshape(Mcrisis,[countries*years-countries*(nlags-1),1]);

constant= ones(length(Ly),1);
tt=1:1:length(Ly)/countries;
for i=1:countries
    Mt(:,i)=tt;
end

trend=reshape(Mt,[countries*years-countries*(nlags-1),1]);
tsq=trend.^2;

DC=zeros((years-nlags+1)*countries,countries);
for i=1:countries
    DC((i-1)*(years-nlags+1)+1:i*(years-nlags+1),i)=1;
end

DT = zeros((years-nlags+1)*countries,years-nlags+1);
for i=1:countries
    DT((i-1)*(years-nlags+1)+1,1)=1;
end

for j=1:years-nlags+1
      DT(j:end,j)=DT(1:end-j+1,1);
end

%Creating part of the new variables that enter the RHS, which are in
%levels, not logs.
Mynew=reshape(y,[years,countries]);
Mynew=Mynew(nlags:end,:);
ynew = reshape(Mynew,[countries*years-countries*(nlags-1),1]);

MLynew=reshape(y,[years,countries]);
MLynew=MLynew(nlags-1:end-1,:);
Lynew = reshape(MLynew,[countries*years-countries*(nlags-1),1]);

Mgnew=reshape(g,[years,countries]);
Mgnew=Mgnew(nlags:end,:);
gnew = reshape(Mgnew,[countries*years-countries*(nlags-1),1]);

MLgnew=reshape(g,[years,countries]);
MLgnew=MLgnew(nlags-1:end-1,:);
Lgnew = reshape(MLgnew,[countries*years-countries*(nlags-1),1]);


%Gshock identification
egy=0.17;
ge=lg-(egy*ly);
Mshock=reshape(ge,[years,countries]);
Mshock=Mshock(nlags:end,:);
shock = reshape(Mshock,[countries*years-countries*(nlags-1),1]);

hor=5;
zz=zeros(1,hor);

I_a = (1-Iexregime).*Ir;
I_b=(1-Iexregime).*(1-Ir);
I_c=Iexregime;

MIa = reshape(I_a,[years,countries]);
MIa(1:nlags-1,:) =[];
Ia= reshape(MIa,[countries*years-countries*(nlags-1),1]);

MIb = reshape(I_b,[years,countries]);
MIb(1:nlags-1,:) =[];
Ib= reshape(MIb,[countries*years-countries*(nlags-1),1]);

MIc = reshape(I_c,[years,countries]);
MIc(1:nlags-1,:) =[];
Ic= reshape(MIc,[countries*years-countries*(nlags-1),1]);

%Vector of conrols
X=[Ia.*shock, Ib.*shock, Ic.*shock, constant, Ia, Ib,... 
  Ia.*Lg, Ia.*L2g, Ia.*Ly, Ia.*L2y,...
   Ia.*Lr, Ia.*L2r, Ia.*Linfl, Ia.*L2infl, Ia.*Lex, Ia.*L2ex, Ia.*Iwar,... 
     Ib.*Lg, Ib.*L2g, Ib.*Ly, Ib.*L2y, ...
     Ib.*Lr, Ib.*L2r, Ib.*Linfl, Ib.*L2infl, Ib.*Lex, Ib.*L2ex, Ib.*Iwar,...
      Ic.*Lg, Ic.*L2g, Ic.*Ly, Ic.*L2y,... 
      Ic.*Lr, Ic.*L2r, Ic.*Linfl, Ic.*L2infl, Ic.*Lex, Ic.*L2ex,...
    DT(:,1:end-hor), DC(:,2:end)];
rpos=1; %position of the shock
[r,nn]=size(X);

t=size(X,1)/countries;
Mc = reshape(country,[t,countries]);


%Vector of conrols
X=[Ia.*shock, Ib.*shock, Ic.*shock, constant, Ia, Ib,... 
  Ia.*Lg, Ia.*L2g, Ia.*Ly, Ia.*L2y,...
   Ia.*Lr, Ia.*L2r, Ia.*Linfl, Ia.*L2infl, Ia.*Lex, Ia.*L2ex, Ia.*Iwar,... 
     Ib.*Lg, Ib.*L2g, Ib.*Ly, Ib.*L2y,...
     Ib.*Lr, Ib.*L2r, Ib.*Linfl, Ib.*L2infl, Ib.*Lex, Ib.*L2ex, Ib.*Iwar,...
      Ic.*Lg, Ic.*L2g, Ic.*Ly, Ic.*L2y,... 
      Ic.*Lr, Ic.*L2r, Ic.*Linfl, Ic.*L2infl, Ic.*Lex, Ic.*L2ex,...
    DT(:,1:end-hor), DC(:,2:end)];
rpos=1; %position of the shock
[r,nn]=size(X);

t=size(X,1)/countries;
Mc = reshape(country,[t,countries]);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%LINEAR (Point-to-Point)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear reg b bint se
clear confidence

ggnc=zeros(size(gnew,1),hor);
gIVc=zeros(size(ggnc,1),hor);
ggnd=zeros(size(gnew,1),hor);
gIVd=zeros(size(ggnd,1),hor);

xc=[Ia.*shock, Ib.*shock, Ic.*shock, constant, Ia, Ib,... 
  Ia.*Lg, Ia.*L2g, Ia.*Ly, Ia.*L2y,...
   Ia.*Lr, Ia.*L2r, Ia.*Linfl, Ia.*L2infl, Ia.*Lex, Ia.*L2ex, Ia.*Iwar,... 
     Ib.*Lg, Ib.*L2g, Ib.*Ly, Ib.*L2y,...
     Ib.*Lr, Ib.*L2r, Ib.*Linfl, Ib.*L2infl, Ib.*Lex, Ib.*L2ex, Ib.*Iwar,...
      Ic.*Lg, Ic.*L2g, Ic.*Ly, Ic.*L2y,... 
      Ic.*Lr, Ic.*L2r, Ic.*Linfl, Ic.*L2infl, Ic.*Lex, Ic.*L2ex,...
    DT(:,1:end-hor), DC(:,2:end)];

Xc=[constant, Ia, Ib,... 
    Ia.*Lg, Ia.*L2g, Ia.*Ly, Ia.*L2y,...
   Ia.*Lr, Ia.*L2r, Ia.*Linfl, Ia.*L2infl, Ia.*Lex, Ia.*L2ex, Ia.*Iwar,... 
     Ib.*Lg, Ib.*L2g, Ib.*Ly, Ib.*L2y,...
     Ib.*Lr, Ib.*L2r, Ib.*Linfl, Ib.*L2infl, Ib.*Lex, Ib.*L2ex, Ib.*Iwar,...
      Ic.*Lg, Ic.*L2g, Ic.*Ly, Ic.*L2y,... 
      Ic.*Lr, Ic.*L2r, Ic.*Linfl, Ic.*L2infl, Ic.*Lex, Ic.*L2ex,...
    DT(:,1:end-hor), DC(:,2:end)];
rpos=1;


for i=1:hor

y11=(Mgnew(i:end,:) - MLgnew(1:end-i+1,:))./ MLynew(1:end-i+1,:);
y1 = reshape(y11,[countries*(years-(nlags-1)-i+1),1]);
Mcc = Mc(i:end,:);
c0 = reshape(Mcc,[(t-i+1)*countries,1]);

for k=1:countries 

Mx = zeros((t-i+1)*countries,nn);    
Mxx = reshape(X((k-1)*t+1:k*t,:),t, nn);
Mx = Mxx(1:t-i+1,:);
xx((k-1)*(t-i+1)+1:k*(t-i+1),:) = Mx;
end

results = clusterreg(y1, xx(1:end,:), c0);
bg(:,i) = results(:,1);
gIVc(1:end-i+1,i)=xc(1:end-i+1,:)*bg(:,i);
clear xx
end

clear reg b bint se yhat
clear confidence

shockIVc=gIVc.*Ia;
shockIVd=gIVc.*Ib;
shockIVe=gIVc.*Ic;

for i=1:hor

y22=(Mynew(i:end,:) - MLynew(1:end-i+1,:))./ MLynew(1:end-i+1,:);
y2 = reshape(y22,[countries*(years-(nlags-1)-i+1),1]);
Mcc = Mc(i:end,:);
c0 = reshape(Mcc,[(t-i+1)*countries,1]);
gshockc=shockIVc(:,i);
gshockd=shockIVd(:,i);
gshocke=shockIVe(:,i);

for k=1:countries 

xYstatec=[gshockc, gshockd, gshocke, Xc];    
Mx = zeros((t-i+1)*countries,nn);    
Mxx = reshape(xYstatec((k-1)*t+1:k*t,:),t, nn);
Mx = Mxx(1:t-i+1,:);
xx((k-1)*(t-i+1)+1:k*(t-i+1),:) = Mx;
end

results = clusterreg(y2, xx(1:end,:), c0);
by(:,i) = results(:,1);
bint(:,1)= by(:,i)-(results(:,2)*ci);
bint(:,2) = by(:,i)+(results(:,2)*ci);
bint1(:,1)= by(:,i)-(results(:,2)*ci1);
bint1(:,2) = by(:,i)+(results(:,2)*ci1);
confidencey(:,:,i)=bint;
confidencey1(:,:,i)=bint1;

clear xx
end
confidencey=reshape(confidencey,2*nn,hor);
confidencey1=reshape(confidencey1,2*nn,hor);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%NON-LINEAR (Cumulative)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ee=zeros(size(X,1),hor);

ggnc=zeros(size(gnew,1),hor);
yync=zeros(size(gnew,1),hor);
gIVccum=zeros(size(ggnc,1),hor);
ggnd=zeros(size(gnew,1),hor);
yynd=zeros(size(gnew,1),hor);
gIVdcum=zeros(size(ggnc,1),hor);


xa=[Ia.*shock, Ib.*shock, Ic.*shock, Xc];

Mggnc=zeros(years-nlags+1,countries,hor);
for j=1:hor
    Mggn=zeros(years-nlags+1,countries);
    Mggn(1:end-j+1,:)=(Mgnew(j:end,:)- MLgnew(1:end-j+1,:))./MLynew(1:end-j+1,:);
    ggnc(1:end,j)=reshape(Mggn,[countries*(years-(nlags-1)),1]);
    Mggnc(1:end-j+1,:,j)=Mggn(1:end-j+1,:);
end

gg1=zeros(size(ggnc,1),hor);
gg1(1:end,1)=ggnc(1:end,1);
Mgg1=zeros(years-nlags+1,countries,hor);
Mgg1(1:end,:,1)=Mggnc(1:end,:,1);
for j=2:hor
    Mgg1(1:end-j+1,:,j)=Mgg1(1:end-j+1,:,j-1)+Mggnc(1:end-j+1,:,j);
    gg1(1:end,j)=reshape(Mgg1(:,:,j),[countries*(years-(nlags-1)),1]);
end

rpos=1;

for i=1:hor
Myg=reshape(gg1(:,i),[years-nlags+1,countries]);
y33=Myg(1:end-i+1,:);
y3 = reshape(y33,[countries*(years-(nlags-1)-i+1),1]);
Mcc = Mc(i:end,:);
c0 = reshape(Mcc,[(t-i+1)*countries,1]);

for k=1:countries 

Mx = zeros((t-i+1)*countries,nn);    
Mxx = reshape(xa((k-1)*t+1:k*t,:),t, nn);
Mx = Mxx(1:t-i+1,:);
xx((k-1)*(t-i+1)+1:k*(t-i+1),:) = Mx;
end

results = clusterreg(y3, xx(1:end,:), c0);
bg_ccum(:,i) = results(:,1);
se_ccum(:,i) = results(:,2);
f_ccum(:,i) = (results(:,3)).^2;
gIVccum(1:end-i+1,i)=xa(1:end-i+1,:)*bg_ccum(:,i);
clear xx
end

clear reg b bint se yhat
clear confidence


shockIVc=gIVccum.*Ia;
shockIVd=gIVccum.*Ib;
shockIVe=gIVccum.*Ic;

Myync=zeros(years-nlags+1,countries,hor);
for j=1:hor
    Myyn=zeros(years-nlags+1,countries);
    Myyn(1:end-j+1,:)=(Mynew(j:end,:)- MLynew(1:end-j+1,:))./MLynew(1:end-j+1,:);
    yync(1:end,j)=reshape(Myyn,[countries*(years-(nlags-1)),1]);
    Myync(1:end-j+1,:,j)=Myyn(1:end-j+1,:);
end

yy1=zeros(size(yync,1),hor);
yy1(1:end,1)=yync(1:end,1);
Myy1=zeros(years-nlags+1,countries,hor);
Myy1(1:end,:,1)=Myync(1:end,:,1);
for j=2:hor
    Myy1(1:end-j+1,:,j)=Myy1(1:end-j+1,:,j-1)+Myync(1:end-j+1,:,j);
    yy1(1:end,j)=reshape(Myy1(:,:,j),[countries*(years-(nlags-1)),1]);
end

for i=1:hor
Myy=reshape(yy1(:,i),[years-nlags+1,countries]);
y44=Myy(1:end-i+1,:);
y4 = reshape(y44,[countries*(years-(nlags-1)-i+1),1]);
Mcc = Mc(i:end,:);
c0 = reshape(Mcc,[(t-i+1)*countries,1]);
gshockc=shockIVc(:,i);
gshockd=shockIVd(:,i);
gshocke=shockIVe(:,i);


for k=1:countries 

xYstatec=[gshockc, gshockd, gshocke, Xc];    
Mx = zeros((t-i+1)*countries,nn);    
Mxx = reshape(xYstatec((k-1)*t+1:k*t,:),t, nn);
Mx = Mxx(1:t-i+1,:);
xx((k-1)*(t-i+1)+1:k*(t-i+1),:) = Mx;
end

results = clusterreg(y4, xx(1:end,:), c0);
by_cum(:,i) = results(:,1);
se_cum(:,i) = results(:,2);
bint(:,1)= by_cum(:,i)-(results(:,2)*ci);
bint(:,2) = by_cum(:,i)+(results(:,2)*ci);
bint1(:,1)= by_cum(:,i)-(results(:,2)*ci1);
bint1(:,2) = by_cum(:,i)+(results(:,2)*ci1);
confidencey_cum(:,:,i)=bint;
confidencey_cum1(:,:,i)=bint1;

clear xx
end
confidencey_cum=reshape(confidencey_cum,2*nn,hor);
confidencey_cum1=reshape(confidencey_cum1,2*nn,hor);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DIFFERENCE (cumulative)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ee=zeros(size(X,1),hor);

ggnc=zeros(size(gnew,1),hor);
yync=zeros(size(gnew,1),hor);
gIVccum=zeros(size(ggnc,1),hor);
ggnd=zeros(size(gnew,1),hor);
yynd=zeros(size(gnew,1),hor);
gIVdcum=zeros(size(ggnc,1),hor);


xa=[Ia.*shock, Ib.*shock, Ic.*shock, Xc];

Mggnc=zeros(years-nlags+1,countries,hor);
for j=1:hor
    Mggn=zeros(years-nlags+1,countries);
    Mggn(1:end-j+1,:)=(Mgnew(j:end,:)- MLgnew(1:end-j+1,:))./MLynew(1:end-j+1,:);
    ggnc(1:end,j)=reshape(Mggn,[countries*(years-(nlags-1)),1]);
    Mggnc(1:end-j+1,:,j)=Mggn(1:end-j+1,:);
end

gg1=zeros(size(ggnc,1),hor);
gg1(1:end,1)=ggnc(1:end,1);
Mgg1=zeros(years-nlags+1,countries,hor);
Mgg1(1:end,:,1)=Mggnc(1:end,:,1);
for j=2:hor
    Mgg1(1:end-j+1,:,j)=Mgg1(1:end-j+1,:,j-1)+Mggnc(1:end-j+1,:,j);
    gg1(1:end,j)=reshape(Mgg1(:,:,j),[countries*(years-(nlags-1)),1]);
end

rpos=1;

for i=1:hor
Myg=reshape(gg1(:,i),[years-nlags+1,countries]);
y33=Myg(1:end-i+1,:);
y3 = reshape(y33,[countries*(years-(nlags-1)-i+1),1]);
Mcc = Mc(i:end,:);
c0 = reshape(Mcc,[(t-i+1)*countries,1]);

for k=1:countries 

Mx = zeros((t-i+1)*countries,nn);    
Mxx = reshape(xa((k-1)*t+1:k*t,:),t, nn);
Mx = Mxx(1:t-i+1,:);
xx((k-1)*(t-i+1)+1:k*(t-i+1),:) = Mx;
end

results = clusterreg(y3, xx(1:end,:), c0);
bg_ccum(:,i) = results(:,1);
gIVccum(1:end-i+1,i)=xa(1:end-i+1,:)*bg_ccum(:,i);
clear xx
end

clear reg b bint se yhat
clear confidence


shockIVc=gIVccum;
shockIVd=gIVccum.*Ia;

Myync=zeros(years-nlags+1,countries,hor);
for j=1:hor
    Myyn=zeros(years-nlags+1,countries);
    Myyn(1:end-j+1,:)=(Mynew(j:end,:)- MLynew(1:end-j+1,:))./MLynew(1:end-j+1,:);
    yync(1:end,j)=reshape(Myyn,[countries*(years-(nlags-1)),1]);
    Myync(1:end-j+1,:,j)=Myyn(1:end-j+1,:);
end

yy1=zeros(size(yync,1),hor);
yy1(1:end,1)=yync(1:end,1);
Myy1=zeros(years-nlags+1,countries,hor);
Myy1(1:end,:,1)=Myync(1:end,:,1);
for j=2:hor
    Myy1(1:end-j+1,:,j)=Myy1(1:end-j+1,:,j-1)+Myync(1:end-j+1,:,j);
    yy1(1:end,j)=reshape(Myy1(:,:,j),[countries*(years-(nlags-1)),1]);
end

for i=1:hor
Myy=reshape(yy1(:,i),[years-nlags+1,countries]);
y44=Myy(1:end-i+1,:);
y4 = reshape(y44,[countries*(years-(nlags-1)-i+1),1]);
Mcc = Mc(i:end,:);
c0 = reshape(Mcc,[(t-i+1)*countries,1]);
gshockc=shockIVc(:,i);
gshockd=shockIVd(:,i);
gshocke=shockIVe(:,i);

for k=1:countries 

xYstatec=[gshockc, gshockd, gshocke, Xc];    
Mx = zeros((t-i+1)*countries,nn);    
Mxx = reshape(xYstatec((k-1)*t+1:k*t,:),t, nn);
Mx = Mxx(1:t-i+1,:);
xx((k-1)*(t-i+1)+1:k*(t-i+1),:) = Mx;
end

results = clusterreg(y4, xx(1:end,:), c0);
by_dif(:,i) = results(:,1);
se_dif(:,i) = results(:,2);
bint(:,1)= by_dif(:,i)-(results(:,2)*ci);
bint(:,2) = by_dif(:,i)+(results(:,2)*ci);
bint1(:,1)= by_dif(:,i)-(results(:,2)*ci1);
bint1(:,2) = by_dif(:,i)+(results(:,2)*ci1);
confidencey_dif(:,:,i)=bint;
confidencey_dif1(:,:,i)=bint1;

clear xx
end
confidencey_dif=reshape(confidencey_dif,2*nn,hor);
confidencey_dif1=reshape(confidencey_dif1,2*nn,hor);

ZLB=[by_cum(rpos,1) se_cum(rpos,1) by_cum(rpos,3) se_cum(rpos,3) by_cum(rpos,5) se_cum(rpos,5)]'
Outside_ZLB=[by_cum(rpos+1,1) se_cum(rpos+1,1) by_cum(rpos+1,3) se_cum(rpos+1,3) by_cum(rpos+1,5) se_cum(rpos+1,5)]'
Diff=[by_dif(rpos+1,1) se_dif(rpos+1,1) by_dif(rpos+1,3) se_dif(rpos+1,3) by_dif(rpos+1,5) se_dif(rpos+1,5)]'
