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]=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);

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


Mr=reshape(r,[years,countries]);
My=reshape(ly,[years,countries]);
Mg=reshape(lg,[years,countries]);
Mmg=reshape(lmg,[years,countries]);


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

MIr=zeros(size(Mr));


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


Ir=reshape(MIr,[countries*years,1]);
I=Ir;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Calculate Growth Rate
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n=2;

MDy=zeros(size(My));
MDg=zeros(size(Mg));
MDmg=zeros(size(Mmg));

for i=1:countries
    for j=2:size(MDy,1)
    MDy(j,i)=100*(My(j,i)-My(j-1,i));
    MDg(j,i)=100*(Mg(j,i)-Mg(j-1,i));
    MDmg(j,i)=100*(Mmg(j,i)-Mmg(j-1,i));
    end   
end
Dy=reshape(MDy,[countries*years,1]);
DyZLB=Dy.*Ir;
DyZLB(DyZLB==0)=[];
DyZLB(abs(DyZLB)>abs(mean(DyZLB)+n*std(DyZLB)))=[];
DyNoZLB=Dy.*(1-Ir);
DyNoZLB(DyNoZLB==0)=[];
DyNoZLB(abs(DyNoZLB)>abs(mean(DyNoZLB)+n*std(DyNoZLB)))=[];
DyFull=Dy;
DyFull(abs(DyFull)>abs(mean(DyFull)+n*std(DyFull)))=[];

Dg=reshape(MDg,[countries*years,1]);
DgZLB=Dg.*Ir;
DgZLB(DgZLB==0)=[];
DgZLB(abs(DgZLB)>abs(mean(DgZLB)+n*std(DgZLB)))=[];
DgNoZLB=Dg.*(1-Ir);
DgNoZLB(DgNoZLB==0)=[];
DgNoZLB(abs(DgNoZLB)>abs(mean(DgNoZLB)+n*std(DgNoZLB)))=[];
DgFull=Dg;
DgFull(abs(DgFull)>abs(mean(DgFull)+n*std(DgFull)))=[];


Dmg=reshape(MDmg,[countries*years,1]);
Dmg(Dmg>mean(Dmg)+4*std(Dmg));
DmgZLB=Dmg.*Ir;
DmgZLB(DmgZLB==0)=[];
DmgZLB(abs(DmgZLB)>abs(mean(DmgZLB)+n*std(DmgZLB)))=[];
DmgNoZLB=Dmg.*(1-Ir);
DmgNoZLB(DmgNoZLB==0)=[];
DmgNoZLB(abs(DmgNoZLB)>abs(mean(DmgNoZLB)+n*std(DmgNoZLB)))=[];
DmgFull=Dmg;
DmgFull(abs(DmgFull)>abs(mean(DmgFull)+n*std(DmgFull)))=[];

[rFull,pFull]=corrcoef(Dmg,Dg);
[rNoZLB,pNoZLB]=corrcoef(Dmg.*(1-Ir),Dg.*(1-Ir));
[rZLB,pZLB]=corrcoef(Dmg.*Ir,Dg.*Ir);
corrFull=[rFull(2,1) pFull(2,1)];
corrNoZLB=[rNoZLB(2,1) pNoZLB(2,1)];
corrZLB=[rZLB(2,1) pZLB(2,1)];

StatsZLB=[mean(DyZLB) mean(DgZLB) mean(DmgZLB) std(DyZLB) std(DgZLB) std(DmgZLB)];
StatsNoZLB=[mean(DyNoZLB) median(DgNoZLB) mean(DmgNoZLB) std(DyNoZLB) std(DgNoZLB) std(DmgNoZLB)];

stdFull=[std(DmgFull) std(DgFull) std(DyFull)];
stdNoZLB=[std(DmgNoZLB) std(DgNoZLB) std(DyNoZLB)];
stdZLB=[std(DmgZLB) std(DgZLB) std(DyZLB)];


disp('Standard deviations')
disp('    gM          g         y')
disp('Full')
disp(stdFull);
disp('Outside')
disp(stdNoZLB);
disp('ZLB')
disp(stdZLB);
disp('--------------------------')
disp('Correlation gm-g')
disp('Full')
disp(corrFull);
disp('No ZLB')
disp(corrNoZLB);
disp('ZLB')
disp(corrZLB);
disp('--------------------------')
disp('Observations')
disp('Full')
disp(length(I))
disp('No ZLB')
disp(length(find(I==0)))
disp('ZLB')
disp(length(find(I==1)))




