% MENA Small open economy model (3 good, 2 production)
% written and revised by Hideaki Hirata

clear all % Tomorrow is another day

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Preparation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%howmany = 1000;   % For monte carlo simulation
howmany = 10;     % For VD or impulses
n_var = 28; 
n_exo = 7; 
SSigma = eye(n_exo);
Karajan =[]; Kleiber =[]; Solti=[];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Starting sensitivity analyses %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%for abc=0:21   % For all sensitivity analysis shown below
%for abc=[0 22 23 24 25 26 27 28 29 30]    % For selected sensitivity analysis shown below
for abc=0;      % Just for benchmark

if abc==0       % Benchmark economy
    dp
elseif abc==1      % If productivity shocks drive the economy
    dp
    ssigmar=  1e-10;
    ssigmag=  1e-10; 
    ssigmapx= 1e-10;
    rhor = 0; 
    rhog = 0;
    rhopx= 0; 
elseif abc==2  % If tot shocks drive the economy
    dp
    ssigmaax= 1e-10;
    ssigmaan= 1e-10;
    ssigmar=  1e-10;
    ssigmag=  1e-10; 
    rhox = 0;
    rhon = 0;
    rhor = 0; 
    rhog = 0;
elseif abc==3  % If volatility of productivity shocks increase by 20% 
    dp
    ssigmaax= 1.2*ssigmaax;
    ssigmaan= 1.2*ssigmaan;
elseif abc==4  % If persistence of productivity shocks increase by 20% 
    dp
    rhox= 1.2*rhox;
    rhon= 1.2*rhon;
elseif abc==5  % If volatility of price shocks increase by 20% 
    dp
    ssigmapx= 1.2*ssigmapx;
elseif abc==6  % If persistence of price shocks increase by 20% 
    dp
    rhopx= 1.2*rhopx;
elseif abc==7  % If volatility of interest shocks increase by 20% 
    dp
    ssigmar = 1.2*ssigmar;
elseif abc==8  % If persistence of interest shocks increase by 20% 
    dp
    rhor = 1.2*rhor;
elseif abc==9  % If elastisity of capital adjustment cost in x increase by 20% 
    dp
    etax=etax*1.2;
elseif abc==10 % If elastisity of capital adjustment cost in n increase by 20% 
    dp
    etan=etan*1.2;
elseif abc==11 % If trade balance is same as Israel (-1.02%) 
    dp
    nyx=-0.26/-0.0906*(-0.0102);    
elseif abc==12 % If trade balance is same as Jordan (-33.70%) 
    dp
    nyx=-0.26/-0.0906*(-0.3370);    
elseif abc==13 % If risk aversion parameter increases by 20% 
    dp
    sigma=1.2*sigma;
elseif abc==14 % If elasticity of substitution in consumption function is Ostry and Reinhart (1992) 
    dp
    gamma=1/1.07;    
elseif abc==15 % If elasticity of substitution in consumption function is >1 
    dp
    gamma=1/0.5;    
elseif abc==16 % If bond adjustment cost is very small 
    dp
    pphi=0.001;
elseif abc==17 % If export good increases in consumption by 20%
    dp
    bx=1.2*bx;
elseif abc==18 % If import good increases in consumption by 20%
    dp
    bm=1.2*bm;
elseif abc==19 % If elasticity b/w capital and imported intermediate good is unity
    dp
    z=1.000001;
elseif abc==20 % If elasticity b/w capital and imported intermediate good is two
    dp
    z=2;
elseif abc==21 % If no land usage
    dp
    lyn=0.00001; 
elseif abc==22 % big nx 33.7
    dp
    nyx=-0.26/-0.0906*(-0.8);    
elseif abc==23 % big nx 33.7
    dp
    nyx=-0.26/-0.0906*(-0.6);    
elseif abc==24 % big nx 33.7
    dp
    nyx=-0.26/-0.0906*(-0.4);    
elseif abc==25 % big nx 33.7
    dp
    nyx=-0.26/-0.0906*(-0.2);    
elseif abc==26 % big nx 33.7
    dp
    nyx=-0.26/-0.0906*(-0.0);    
elseif abc==27 % big nx 33.7
    dp
    nyx=-0.26/-0.0906*(0.2);    
elseif abc==28 % big nx 33.7
    dp
    nyx=-0.26/-0.0906*(0.4);    
elseif abc==29 % big nx 33.7
    dp
    nyx=-0.26/-0.0906*(0.6);    
elseif abc==30 % big nx 33.7
    dp
    nyx=-0.26/-0.0906*(0.8);    
end    



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% endogenously derived steady state values (ratios) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

byx=-(nyx/px+fyx)/r; % initial asset holding position as percent of pxyx
kyn=beta*alpha3/(1-beta*(1-deltan)); % production n
hyn=lyn^(-alpha2/alpha1)*kyn^(-alpha3/alpha1);
iyn=kyn*deltan;
cyn=1-gyn;

kkk1=omega*(1-z)/z/(1-omega);
aa1=a*beta*(1-omega)/(1-beta*(1-deltax));
aa2=px*(1-a)*(1-omega);
pn=px*omega/alpha1*hyn*(a*aa1^((1-z)/z)+(1-a)*aa2^((1-z)/z))^(1/kkk1);  

hyx=hyn*px*omega/alpha1/pn;
kyx=aa1^(1/z)*hyx^kkk1;
syx=aa2^(1/z)*hyx^kkk1;
iyx=kyx*deltax;

pn1=bn^(1/gamma)*pn^((gamma-1)/gamma)+bx^(1/gamma)*px^((gamma-1)/gamma);
pn2=bm^(1/gamma);
p=(pn1+pn2)^(gamma/(gamma-1));

cxc=bx^(1/gamma)/(px/p)^(1/gamma);
cmc=bm^(1/gamma)/(1/p)^(1/gamma);
cnc=bn^(1/gamma)/(pn/p)^(1/gamma);

cyx1=1-iyx/px-syx/px+fyx+byx*r;
cyx2=1+cmc/cxc/px;
cyx3=1/px*iyn*(bn/bx)^(1/gamma)*(px/pn)^(1/gamma)/cyn;
cyx=cyx1/(cyx2+cyx3);
inyx=cyx3*cyx;

h=(hyx/(omega*bx*cxc^(-gamma)))^(1/(1-theta));

% now ss for all individual variables
%l=(1+hxl+hnl)^(-1);
%hx=hxl*l
%hn=hnl*l

%%%For finding the SS hx and hm%%%%
%for dd=0:0.1:1;
%for dd=0.300:0.001:0.400;
%for dd=0.5430:0.00001:0.5440;
%%%For finding the SS hx and hm%%%%

%%%Once the SS hx and hm are found%%%%
   dd=0.37660;
%%%Once the SS hx and hm are found%%%%

hx=dd*h;
    hn=h-hx;
    ll=1-h;
yn=hn/hyn;
yx=hx/hyx;
y=px*yx+pn*yn;
cn=cyn*yn;
cx=cyx*yx;
cm=cmc/cxc*cx;
cc=px*cx+cm+pn*cn;   % pt*ct
ccc=cc/p;   %ct
in=iyn*yn;
ix=iyx*yx;
i=ix+in; %s is an imported intermediate good, not imported capital good

gn=gyn*yn;
f=fyx*yx;
g=pn*gn;
s=syx*yx;

%Once the SS nx is found
nx=nyx*yx;
%Once the SS nx is found

%%%For finding the SS nx%%%%
%nx=y-cc-g-i-s;
%[dd nx]
%end
%%%For finding the SS nx%%%%

kn=kyn*yn;
kx=kyx*yx;
k=kx+kn;

consx=cx/cc;
consm=cm/cc;
consn=pn*cn/cc;

bb=byx*px*yx;
ca=nx+f+r*bb;
sav=y-cc-g-s;

import=cm+s+in+ix;
export=nx+import;

exportshare=export/y;

%define variable xx
xx=ccc-(1-ll)^theta/theta;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% program for running gensys %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% note that r, B and K at time (t+1) are known at t+1. 
% See footnote 15.
% defining matrices

g0 = zeros(n_var,n_var);  % x(t+1)
g1 = zeros(n_var,n_var);  % x(t)
c = zeros(n_var,1);       % constant
psi = zeros(n_var,n_exo); % shocks(t+1)
pi = zeros(n_var,n_var);  % expectations(t+1)

%variables 

%U

g0(1,1)=1;
g0(1,2)=-(1-sigma)*ccc/xx;
g0(1,10)=-(1-sigma)*ll/xx*(1-ll)^(theta-1);

%c

g0(2,2)=-1;
g0(2,3)=bx*cxc^(1-gamma);
g0(2,4)=bm*cmc^(1-gamma);
g0(2,5)=bn*cnc^(1-gamma);

%l

g0(3,10)=1;
g0(3,11)=hx/ll;
g0(3,12)=hn/ll;


%budget constraint for n

g0(4,7)=1;
g0(4,5)=-cyn;
g0(4,26)=-gyn;

%budget constraint for x

g0(5,6)=1;
g0(5,3)=-cyx;
g0(5,4)=-cyx*(cmc/cxc)/px;;
g0(5,8)=-iyx/px;
g0(5,13)=-syx/px;
g0(5,25)=1;
g1(5,24)=-byx;  % B at t+1 are known at time t+1
g0(5,27)=1-cyx-iyx;
g0(5,16)=-beta*(1+r);
g0(5,9)=-inyx;
g1(5,16)=-(1+r);

% c-cx

g0(6,1)=sigma/(1-sigma);
g0(6,19)=1;
g0(6,27)=1;
g0(6,2)=-gamma;
g0(6,3)=gamma;

%cm-modified

g0(7,27)=1;
g0(7,3)=gamma;
g0(7,4)=-gamma;

%cn-modified

g0(8,17)=1;
g0(8,5)=gamma;
g0(8,4)=-gamma;

% cx

g0(9,2)=1;
g0(9,3)=-1;
g0(9,18)=1/gamma;
g0(9,27)=-1/gamma;

%hx

g0(10,1)=sigma/(1-sigma);
g0(10,19)=1;
g0(10,27)=1;
g0(10,6)=1;
g0(10,11)=-1;
g0(10,10)=(theta-1)*ll/(1-ll);

%hn

g0(11,1)=sigma/(1-sigma);
g0(11,19)=1;
g0(11,17)=1;
g0(11,7)=1;
g0(11,12)=-1;
g0(11,10)=(theta-1)*ll/(1-ll);;

%bond

g0(12,19)=1;
g1(12,19)=1;
g1(12,16)=(pphi*px*yx)/beta/(1+r);
%g1(12,16)=0;
g1(12,24)=-1/(1+r);
pi(12,12)=1;
%pi(12,:)=g0(12,:);

%yx  

g0(13,6)=-1;
g0(13,22)=1;
g0(13,11)=omega;
g1(13,14)=-kyx*(1-beta*(1-deltax))/beta;
g0(13,13)=syx/px;

%yn

g0(14,7)=-1;
g0(14,23)=1;
g0(14,12)=alpha1;
g0(14,28)=alpha2;
g1(14,15)=-alpha3;

%capital accumulation x

g0(15,14)=1;
g1(15,14)=1-deltax;
g0(15,8)=-deltax;

%capital accumulation n

g0(16,15)=1;
g1(16,15)=1-deltan;
g0(16,9)=-deltan;

%ix

g0(17,8)=1;
g1(17,14)=1;
g0(17,20)=-etax;
g0(17,19)=etax;
%g0(17,27)=etax;

%s

g0(18,6)=1;
g0(18,27)=1;
g1(18,14)=kyx*(1-z)*(1-beta*(1-deltax))/beta/(1-omega);
g0(18,13)=-(1-kyx*(1-z)*(1-beta*(1-deltax))/beta/(1-omega));

%kx

g0(19,20)=beta;
g0(19,19)=1-beta;
g0(19,27)=1-beta;
g0(19,6)=1-beta*(1-deltax);
g0(19,13)=-(1-beta*(1-deltax))*((1-z)/(1-omega)*syx/px);
g1(19,14)=(1-beta*(1-deltax))*(1-(1-z)/(1-omega)*syx/px);
g1(19,20)=1;
pi(19,19)=1;
%pi(19,:)=g0(19,:);

%in

g0(20,9)=1;
g1(20,15)=1;
g0(20,21)=-etan;
g0(20,19)=etan;
%g0(20,17)=etan;

%kn

g0(21,21)=beta;
g0(21,19)=1-beta;
g0(21,17)=1-beta;
g0(21,7)=1-beta*(1-deltan);
g1(21,15)=(1-beta*(1-deltan));
g1(21,21)=1;
pi(21,21)=1;
%pi(21,:)=g0(21,:);

%Ax

g0(22,22)=1;
g1(22,22)=rhox;
psi(22,1)=1;

%An

g0(23,23)=1;
g1(23,23)=rhon;
psi(23,2)=1;

%r

g0(24,24)=1;
g1(24,24)=rhor;
psi(24,3)=1;

%f

g0(25,25)=1;
g1(25,25)=rhof;
psi(25,4)=1;

%Gn

g0(26,26)=1;
g1(26,26)=rhog;
psi(26,5)=1;

%px

g0(27,27)=1;
g1(27,27)=rhopx;
psi(27,6)=1;

%Ln

g0(28,28)=1;
g1(28,28)=rholn;
psi(28,7)=1;

[G1,C,impact,fmat,fwt,ywt,gev,eu]=gensys(g0,g1,c,psi,pi);


SSigma(1,1) = ssigmaax^2;
SSigma(2,2) = ssigmaan^2;
SSigma(3,3) = ssigmar^2;
SSigma(4,4) = ssigmaf^2;
SSigma(5,5) = ssigmag^2;
SSigma(6,6) = ssigmapx^2;
SSigma(7,7) = ssigmal^2;

SSigma(1,2) = coraxan*ssigmaax*ssigmaan;
SSigma(1,3) = coraxr*ssigmaax*ssigmar;
SSigma(1,5) = coraxg*ssigmaax*ssigmag;
SSigma(1,6) = coraxpx*ssigmaax*ssigmapx;
SSigma(2,3) = coranr*ssigmaan*ssigmar;
SSigma(2,5) = corang*ssigmaan*ssigmag;
SSigma(2,6) = coranpx*ssigmaan*ssigmapx;
SSigma(3,5) = corrg*ssigmar*ssigmag;
SSigma(3,6) = corrpx*ssigmar*ssigmapx;
SSigma(5,6) = corgpx*ssigmag*ssigmapx;

upper = chol(SSigma); lower = upper';
varsh = diag(lower.^2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

howlong = 41; %1960-2000
period=howlong;
randn('state',123)
vec_moments = [];
persi=[];

%28 var and 7 shocks

for simul = 1:howmany
	[simul howmany];
    simul;
    seri = [];
	added = zeros(n_var,1);
	for pe = 1:howlong
      added = C+G1*added+impact*lower*[randn(1,n_exo)]';
		seri = [seri added];
	end
   
   cxhat = seri(3,:);
   cmhat = seri(4,:);
   cnhat = seri(5,:);
   yxhat = seri(6,:);
   ynhat = seri(7,:);
   ixhat = seri(8,:);
   inhat = seri(9,:);
   hxhat = seri(11,:);
   hnhat = seri(12,:);
   shat = seri(13,:);
   kxhat = seri(14,:);
   knhat = seri(15,:);
   bhat = seri(16,:);
   pnhat = seri(17,:);
   phat = seri(18,:);   
   
   axhat = seri(22,:);
   anhat = seri(23,:);
   rhat = seri(24,:);
   fhat = seri(25,:);
   gnhat = seri(26,:);
   pxhat = seri(27,:);   
   lnhat = seri(28,:);

   %aggregate variable in terms of import good prices: c,y,i,tb,ca
   
   yhat=(px*yx/y)*(pxhat+yxhat)+(pn*yn/y)*(pnhat+ynhat);
    gnpbar=y+r*bb;
    gnphat=y/gnpbar*yhat+px*yx*r/y*bhat+bb/gnpbar*rhat;
   cchat=(px*cx/cc)*(pxhat+cxhat)+(cm/cc)*cmhat+(pn*cn/cc)*(pnhat+cnhat); %p*c, aggregate nominal consumption
   ihat=(ix/i)*ixhat+(in/i)*inhat;
   savhat=y/sav*yhat-cc/sav*cchat-g/sav*(pnhat+gnhat)-(s/sav)*shat;
   intpayhat=px*yx*r/y*bhat+bb/y*rhat; %as percentage of y
   nxhat=(((px*yx)/y)*(pxhat+yxhat))-(((px*cx)/y)*(pxhat+cxhat))-((cm/y)*cmhat+(ix/y)*ixhat+(in/y)*inhat+(s/y)*shat);
   cahat=nxhat(2:period)+px*yx/y*fhat(2:period)+bb/y*rhat(1:period-1)+px*yx*r/y*bhat(1:period-1);  %percentage from y
   cahat=[0,cahat];
   hhat=hx/(hx+hn)*hxhat+ hn/(hx+hn)*hnhat;    
   openhat=(((px*yx)/y)*(pxhat+yxhat))-(((px*cx)/y)*(pxhat+cxhat))+((cm/y)*cmhat+(ix/y)*ixhat+(in/y)*inhat+(s/y)*shat);

   
   %sectoral varaibles
   
   pxyxhat=pxhat+yxhat;  %px*yx
   pxcxhat=pxhat+cxhat;  %px*cx
   pxixhat=ixhat;
   
   pnynhat=pnhat+ynhat; %pn*yn
   pncnhat=pnhat+cnhat; %pn*cn
   pninhat=inhat;

   
   %find persistence from AR(1)
    
    perh2  =[yhat', gnphat', cchat', ihat', hhat',nxhat',savhat', intpayhat', bhat',phat', pxyxhat', pnynhat', openhat']; 
    persis=[];
   for zz=1:size(perh2,2)
       per=[ones(howlong-1,1) perh2(1:howlong-1,zz)]\perh2(2:howlong,zz);
       persis=[persis per(2)];
   end

%H-P filtering: except TB/Y, B/Y, intpay/Y.
    paramm=100;
	cxhat = cxhat-(hpfilter(cxhat,paramm))';
  	cmhat = cmhat-(hpfilter(cmhat,paramm))';
   	cnhat = cnhat-(hpfilter(cnhat,paramm))';
   	yxhat = yxhat-(hpfilter(yxhat,paramm))';
   	ynhat = ynhat-(hpfilter(ynhat,paramm))';
   	ixhat = ixhat-(hpfilter(ixhat,paramm))';
    inhat = inhat-(hpfilter(inhat,paramm))';
    hxhat = hxhat-(hpfilter(hxhat,paramm))';
    hnhat = hnhat-(hpfilter(hnhat,paramm))';
    shat = shat-(hpfilter(shat,paramm))';
    kxhat = kxhat-(hpfilter(kxhat,paramm))';
    knhat = knhat-(hpfilter(knhat,paramm))';
    pnhat = pnhat-(hpfilter(pnhat,paramm))'; %price of nontraded good
    phat = phat-(hpfilter(phat,paramm))';    %real exchange rate
    yhat = yhat-(hpfilter(yhat,paramm))';
    cchat = cchat-(hpfilter(cchat,paramm))';
    ihat = ihat-(hpfilter(ihat,paramm))';
    savhat = savhat-(hpfilter(savhat,paramm))';
    gnphat = gnphat-(hpfilter(gnphat,paramm))';
    pxyxhat = pxyxhat-(hpfilter(pxyxhat,paramm))';   
    pnynhat = pnynhat-(hpfilter(pnynhat,paramm))';   
      
%sectoral variables   
seriix=[yxhat',cxhat',ixhat',hxhat',shat',kxhat']';
seriin=[ynhat',cnhat',inhat',hnhat',knhat']';
seriic=[cxhat',cmhat',cnhat']'; %or seriic=[pxcxhat',cmhat',pncnhat']';

%aggregate varibles 
seriiagg=[yhat',cchat',ihat',nxhat',cahat' ,openhat']';
seriikose=[yhat', gnphat', cchat', ihat', hhat',nxhat',savhat', intpayhat', bhat',phat',openhat']'; %including new variables
seriih2  =[yhat', gnphat', cchat', ihat', hhat',nxhat',savhat', intpayhat', bhat',phat', pxyxhat', pnynhat',openhat']'; %including new variables
    
%price variables
seriip=[pxhat',pnhat',phat']';
seriiyp=[yhat',cchat',nxhat',cahat',pxhat',pnhat',phat' ,openhat']';

%all variables--y,cc,i,nx,ca,px,pn,p,yx,cx,ix,hx,s,kx,yn,cn,in,hn,kn,cm (20*20 matrix)
seriiall=[seriiagg',seriip',seriix',seriin',cmhat']';
   
% s.d. and corr matrix of variables

%serii=seriikose;         
serii=seriih2;         

[rownum,colnum]=size(serii);

dum=tril(corrcoef(serii'),-1)+diag(sqrt(diag(cov(serii'))));
vec_moments = [vec_moments vech(dum)];
persi=[persi; persis];   

end

mean_moments = vech_1((mean(vec_moments'))',rownum);

%standard deviation
sd_gdp=mean_moments(1,1);
sd_gnp=mean_moments(2,2);
sd_con=mean_moments(3,3);
sd_inv=mean_moments(4,4);
sd_hour=mean_moments(5,5);
sd_nx=mean_moments(6,6);
sd_sav=mean_moments(7,7);
sd_intpay=mean_moments(8,8);
sd_nfa=mean_moments(9,9);
sd_rer=mean_moments(10,10);
sd_pxyx=mean_moments(11,11);
sd_pnyn=mean_moments(12,12);
sd_open=mean_moments(13,13);

%relative standard deviation
rsd_gdp=mean_moments(1,1)/mean_moments(1,1);
rsd_gnp=mean_moments(2,2)/mean_moments(1,1);
rsd_c=mean_moments(3,3)/mean_moments(1,1); 
rsd_i=mean_moments(4,4)/mean_moments(1,1);
rsd_h=mean_moments(5,5)/mean_moments(1,1);
rsd_nx=mean_moments(6,6)/mean_moments(1,1);
rsd_sav=mean_moments(7,7)/mean_moments(1,1);
rsd_interpay=mean_moments(8,8)/mean_moments(1,1);
rsd_nfa=mean_moments(9,9)/mean_moments(1,1);
rsd_rer=mean_moments(10,10)/mean_moments(1,1);
rsd_pxyx=mean_moments(11,11)/mean_moments(1,1);
rsd_pnyn=mean_moments(12,12)/mean_moments(1,1);
rsd_open=mean_moments(13,13)/mean_moments(1,1);

%correlation with gdp
cor_y_gdp=mean_moments(1,1);
cor_y_gnp=mean_moments(1,2);
cor_y_c=mean_moments(1,3);
cor_y_i=mean_moments(1,4);
cor_y_h=mean_moments(1,5);
cor_y_nx=mean_moments(1,6);
cor_y_sav=mean_moments(1,7);
cor_y_interpay=mean_moments(1,8);
cor_y_nfa=mean_moments(1,9);
cor_y_rer=mean_moments(1,10);
cor_y_pxyx=mean_moments(1,11);
cor_y_pnyn=mean_moments(1,12);
cor_y_open=mean_moments(1,13);
cor_s_i=mean_moments(4,7);

%persistence
persistence=mean(persi);
   per_gdp=persistence(1);
   per_gnp=persistence(2);
   per_con=persistence(3);
   per_inv=persistence(4);
   per_hour=persistence(5);
   per_nx=persistence(6);
   per_sav=persistence(7);
   per_interpay=persistence(8);
   per_nfa=persistence(9);
   per_rer=persistence(10);
   per_pxyx=persistence(11);
   per_pnyn=persistence(12);
   per_open=persistence(13);

AAA=[sd_gdp	    rsd_gdp	    1	           per_gdp	
     sd_pxyx	rsd_pxyx	cor_y_pxyx	   per_pxyx	
     sd_pnyn	rsd_pnyn	cor_y_pnyn	   per_pnyn	
     sd_con	    rsd_c       cor_y_c	       per_con	
     sd_inv	    rsd_i	    cor_y_i	       per_inv	
	 sd_nx      rsd_nx      cor_y_nx       per_nx	
	 sd_open    rsd_open    cor_y_open     per_open];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Variance decomposition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

vvc_cx 	=[]; vvc_cm  =[]; vvc_cn    =[]; vvc_yx 	=[]; vvc_yn =[];
vvc_ix 	=[]; vvc_in  =[]; vvc_hx    =[]; vvc_hn 	=[]; vvc_b 	=[];
vvc_y 	=[]; vvc_c   =[]; vvc_i     =[]; vvc_nx 	=[]; vvc_h 	=[];
vvc_open=[];

aa=zeros(6,n_var);

%Y
aa(1,7)=pn*yn/y;
aa(1,6)=px*yx/y;
aa(1,17)=pn*yn/y;
aa(1,27)=px*yx/y;

%C
aa(2,3)=px*cx/cc;
aa(2,4)=cm/cc;
aa(2,5)=pn*cn/cc;
aa(2,17)=pn*cn/cc;
aa(2,27)=px*cx/cc;

%i
aa(3,8)=ix/i;
aa(3,9)=in/i;

aa(4,27)=px*yx/y-px*cx/y; % pxhat
aa(4,06)=px*yx/y;         % yxhat
aa(4,03)=-px*cx/y;        % cxhat
aa(4,04)=-cm/y;    % cmhat
aa(4,08)=-ix/y;    % ixhat
aa(4,09)=-in/y;    % inhat
aa(4,13)=-s/y ;    % shat

%h
aa(5,11)=hx/h;
aa(5,12)=hx/h;

%open
aa(6,27)=px*yx/y-px*cx/y; % pxhat
aa(6,06)=px*yx/y;         % yxhat
aa(6,03)=-px*cx/y;        % cxhat
aa(6,04)=cm/y;    % cmhat
aa(6,08)=ix/y;    % ixhat
aa(6,09)=in/y;    % inhat
aa(6,13)=s/y ;    % shat




for mmm=1:10
finalpe=mmm;

%for individual variables
mnum=size(impact,2);
xx=zeros(n_var,mnum);

for iiii=0:finalpe-1
    xxx = (G1^iiii*impact).^2; 
    xx=xx+xxx;
end

totalvar=xx*varsh;

var1=xx(:,1)*varsh(1);
var2=xx(:,2)*varsh(2);
var3=xx(:,3)*varsh(3);
var4=xx(:,4)*varsh(4);
var5=xx(:,5)*varsh(5);
var6=xx(:,6)*varsh(6);
var7=xx(:,7)*varsh(7);

varr=[var1 var2 var3 var4 var5 var6 var7];

var_cx=varr(3,:)./totalvar(3);
var_cm=varr(4,:)./totalvar(4);
var_cn=varr(5,:)./totalvar(5);
var_yx=varr(6,:)./totalvar(6);
var_yn=varr(7,:)./totalvar(7);
var_ix=varr(8,:)./totalvar(8);
var_in=varr(9,:)./totalvar(9);
var_hx=varr(11,:)./totalvar(11);
var_hn=varr(12,:)./totalvar(12);
var_b=varr(16,:)./totalvar(16);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%for aggregate variables

zz=zeros(6,mnum);

for iiiiiii=0:finalpe-1
    zzz = (aa*G1^iiiiiii*impact).^2; 
    zz=zz+zzz;
end

varsh1=diag(diag(varsh));
totalvarz=zz*varsh1;

var1z=zz(:,1)*varsh(1);
var2z=zz(:,2)*varsh(2);
var3z=zz(:,3)*varsh(3);
var4z=zz(:,4)*varsh(4);
var5z=zz(:,5)*varsh(5);
var6z=zz(:,6)*varsh(6);
var7z=zz(:,7)*varsh(7);

varrz=[var1z var2z var3z var4z var5z var6z var7z];

var_y=varrz(1,:)./totalvarz(1);
var_c=varrz(2,:)./totalvarz(2);
var_i=varrz(3,:)./totalvarz(3);
var_nx=varrz(4,:)./totalvarz(4);
var_h=varrz(5,:)./totalvarz(5);
var_open=varrz(6,:)./totalvarz(6);


vvc_cx 	=[	vvc_cx 	; var_cx 	];
vvc_cm 	=[	vvc_cm 	; var_cm 	];
vvc_cn 	=[	vvc_cn 	; var_cn 	];
vvc_yx 	=[	vvc_yx 	; var_yx 	];
vvc_yn 	=[	vvc_yn 	; var_yn 	];
vvc_ix 	=[	vvc_ix 	; var_ix 	];
vvc_in 	=[	vvc_in 	; var_in 	];
vvc_hx 	=[	vvc_hx 	; var_hx 	];
vvc_hn 	=[	vvc_hn 	; var_hn 	];
vvc_b 	=[	vvc_b 	; var_b 	];
vvc_y 	=[	vvc_y 	; var_y 	];
vvc_c 	=[	vvc_c 	; var_c 	];
vvc_i 	=[	vvc_i 	; var_i 	];
vvc_nx 	=[	vvc_nx 	; var_nx 	];
vvc_h 	=[	vvc_h 	; var_h 	];
vvc_open=[	vvc_open; var_open 	];

end

vdecall=[	vvc_y 	
	vvc_yx 	
	vvc_yn 	
	vvc_c 	
    vvc_cx 	
	vvc_cm 	
	vvc_cn 	
	vvc_i 	
	vvc_ix 	
	vvc_in 	
	vvc_h
    vvc_hx 	
	vvc_hn 	
	vvc_open 	
	vvc_b 	
    vvc_nx];

uchu=[h nyx byx export+import vvc_y(1,6)];


Solti   =[Solti; uchu];
Karajan =[Karajan; AAA]; 
Kleiber =[Kleiber; vdecall]; 

abc=abc+1;
end;

delete TABLE_sm.txt
delete TABLE_vd.txt
diary ('TABLE_sm.txt')
Karajan
Solti
diary off;
diary ('TABLE_vd.txt')
Kleiber
diary off;