%
% Dynare code for the Devereux-Sutherland solution for asset holdings in an 
% example model.
%
% Written by Alan Sutherland, 23 January 2014
%
% The solution procedure is implemented in two stages. Zero-order portfolio 
% holdings are calculated in Stage 1. First-order portfolio holdings and 
% second and third-order excess returns are calculated in Stage 2.
%

% Variable and parameter declarations

var YK1 YK2 YL1 YL2 M1 M2 NFA Z1 Z2 Z3 Z4 Y1 Y2 C1 C2 P1 P2 cg cd rb1 rb2 re1 re2;

varexo ep1 ep2 ep3 ep4 ep5 ep6 zeta;

parameters BT DK DL DM RH FF;

% Parameter values

BT = 0.98;          
DK = 0.95;
DL = 0.95;
DM = 0.7;
RH = 2.0; 
FF = 0.36;

% The menu of assets and the NFA equation.
% af1, af2, af3 etc are the zero-order portfolio holdings calculated by the
% Devereux-Sutherland (DS) method. af1 etc are set to zero in Stage 1. In
% Stage 2 af1 etc are set to the values calculated in Stage 1. Note that in the 
% DS method zero-order portfolio holdings are defined as alpha/(beta*Y_bar)
% and net foreign assets are defined as NFA/Y_bar, so when the portfolio 
% return and zeta terms are appended to the non-approximated form of the 
% NFA equation af1 etc are % multiplied by beta*Y_bar and zeta is multipled 
% by Y_bar.
% The cell array dr contains a list of return names which correspond to the
% menu of assets.  The final return in the list is used as the numeraire
% asset for defining excess returns.

dr={'re1','re2','rb1','rb2'};                  

parameters af1 af2 af3;

af1 = 0;
af2 = 0;
af3 = 0;

model; 

NFA = exp(rb2)*NFA(-1) + exp(Y1) - exp(C1) + exp(STEADY_STATE(Y1))*( BT*af1*(exp(re1) - exp(rb2)) + BT*af2*(exp(re2) - exp(rb2)) + BT*af3*(exp(rb1) - exp(rb2)) + zeta );

end;

% Main model equations

model;

cd = C1 - C2; 

cg = (1/2)*(C1 + C2);

YK1 = DK*YK1(-1) + ep1;                                                       
YK2 = DK*YK2(-1) + ep2; 
YL1 = DL*YL1(-1) + ep3;                                                       
YL2 = DL*YL2(-1) + ep4;
M1 = DM*M1(-1) + ep5;                                                       
M2 = DM*M2(-1) + ep6;

exp(Y1) = FF*exp(YK1) + (1-FF)*exp(YL1);
exp(Y2) = FF*exp(YK2) + (1-FF)*exp(YL2);
 
exp(Y1) + exp(Y2) = exp(C1) + exp(C2);  

exp(M1)/exp(P1) = exp(Y1);
exp(M2)/exp(P2) = exp(Y2);             

exp(rb1) =  1/(exp(P1)*exp(Z1(-1))); 
exp(rb2) =  1/(exp(P2)*exp(Z2(-1)));

exp(re1) =  exp(YK1)/exp(Z3(-1)); 
exp(re2) =  exp(YK2)/exp(Z4(-1));

BT*( exp(-RH*C2(+1)) * exp(rb2(+1)) ) = exp(-RH*C2);

BT*( exp(-RH*C1(+1)) * exp(rb1(+1)) ) = exp(-RH*C1);
BT*( exp(-RH*C1(+1)) * exp(rb2(+1)) ) = exp(-RH*C1);

BT*( exp(-RH*C1(+1)) * exp(re1(+1)) ) = exp(-RH*C1);
BT*( exp(-RH*C1(+1)) * exp(re2(+1)) ) = exp(-RH*C1);

end;

% Steady state values

steady_state_model;
cd = 0;
cg = 0;
NFA = 0;
Y1 = 0;
Y2 = 0;
YK1 = 0;
YK2 = 0;
YL1 = 0;
YL2 = 0;
M1 = 0;
M2 = 0;
C1 = 0;
C2 = 0;
P1 = 0;
P2 = 0;
rb1 = log(1/BT);
rb2 = log(1/BT);
re1 = log(1/BT);
re2 = log(1/BT);
Z1 = log(BT);
Z2 = log(BT);
Z3 = log(BT);
Z4 = log(BT);
end; 

steady;

% Shock variances

shocks;
var ep1 = 1;
var ep2 = 1;
var ep3 = 1;
var ep4 = 1;
var ep5 = 1;
var ep6 = 1;
end;

% Stage 1: calculate zero-order asset holdings

stoch_simul(order=1,noprint,nomoments,irf=0);

BB=oo_.dr.ghu;
nvar=M_.endo_nbr;
nu=M_.exo_nbr-1;
ordo=oo_.dr.order_var;

for i=1:nvar
    varo{i,1} = M_.endo_names(ordo(i),:);
end

SIGMA=M_.Sigma_e(1:nu,1:nu);

icd=inx('cd',varo);

[dm,na]=size(dr);

for i=1:na
    ia(i)=inx(dr{i},varo);
end

D1=BB(icd,nu+1);
D2=BB(icd,1:nu);

for k=1:na-1
    R1(k,:)=BB(ia(k),nu+1)-BB(ia(na),nu+1);
    R2(k,:)=BB(ia(k),1:nu)-BB(ia(na),1:nu);
end

alpha_tilde=inv(R2*SIGMA*D2'*R1'-R2*SIGMA*R2'*D1)*R2*SIGMA*D2';

% Set af1, af2 etc to calculated values 

npr=M_.param_nbr-na+1;

for i=1:na-1
    alval{i,1}=['M_.params( ' num2str(i+npr) ' )=alpha_tilde(' num2str(i) ',1);'];
    alval{i,2}=['af' num2str(i) '=alpha_tilde(' num2str(i) ',1);'];
end   

for i=1:na-1
    eval(alval{i,1});
    eval(alval{i,2});
end  

% Stage 2: calculate portfolio (first-order) and excess return (third-order) dynamics 
% The gamma vector

stoch_simul(order=2,noprint,nomoments,irf=0);

BB=oo_.dr.ghu;
EE=oo_.dr.ghxu;

nx=M_.npred;

D1=BB(icd,nu+1);
D2=BB(icd,1:nu);

for i=1:nx;
    for j=1:nu
            D5(i,j)=EE(icd,(i-1)*(nu+1)+j);   
    end    
end;

for k=1:na-1
    R2(k,:)=BB(ia(k),1:nu)-BB(ia(na),1:nu);
    R2k=BB(ia(k),1:nu)-BB(ia(na),1:nu);
    for i=1:nx;
        for j=1:nu
                R5k(i,j)=EE(ia(k),(i-1)*(nu+1)+j)-EE(ia(na),(i-1)*(nu+1)+j);
        end    
    end;
    CX2(k,:)=R2k*SIGMA*D5'+D2*SIGMA*R5k'; 
end

VX2=R2*SIGMA*R2'*D1;

gamma=-inv(VX2)*CX2;

% The delta vector

icg=inx('cg',varo);

G1=BB(icg,nu+1);
G2=BB(icg,1:nu);

for i=1:nx;
    for j=1:nu
            G5(i,j)=EE(icg,(i-1)*(nu+1)+j);   
    end    
end;

rs=RH*G2*SIGMA*R2';              % second-order excess returns

Gx5=G5+G1*gamma'*R2;

AA=oo_.dr.ghx;
Rb3=AA(ia(na),:);

for k=1:na-1
    R1k=BB(ia(k),nu+1)-BB(ia(na),nu+1);
    R2k=BB(ia(k),1:nu)-BB(ia(na),1:nu);
    for i=1:nx;
        for j=1:nu
            R5k(i,j)=EE(ia(k),(i-1)*(nu+1)+j)-EE(ia(na),(i-1)*(nu+1)+j);
        end    
    end;
    Rx5k=R5k+R1k*gamma'*R2;
    delta(k,:)=RH*R2k*SIGMA*Gx5'+RH*G2*SIGMA*Rx5k'+rs(k)*Rb3;
end

% First order impulse responses (third order for excess returns)

nsk=M_.nstatic;

AA=[zeros(nvar,nsk) AA zeros(nvar,nvar-nx-nsk)];

hr=50;

for j=1:nu
    
    clear IX;
    
    SVX=zeros(nu+1,1);
    SVX(j,1)=1;
    
    YT=BB*SVX;
    IX(:,1)=[1; YT];
    
    for t=2:hr+1
        YT=AA*YT;
        IX(:,t)=[t; YT];
    end
    IX=IX';
    
    stv=IX(:,nsk+2:nx+nsk+1);
    
    afv=(gamma*stv')';      % asset holdings (first order)
    rxv=(delta*stv')';      % excess returns (third order)
    
    IX=[IX afv rxv];
    
    for i=1:na-1
        varo{nvar+i}=['af' num2str(i)];
    end
    for i=1:na-1
        varo{nvar+na-1+i}=['ex' num2str(i)];
    end
    
    ZZ{1,1}='t';
    for i=1:nvar+2*(na-1)
        ZZ{1,i+1}=varo{i};
    end
    
    for i=1:nvar+1+2*(na-1)
        for t=1:hr
            ZZ{t+1,i}=IX(t,i);
        end
    end
    
    IPR{j,1}=ZZ;

end

% Impulse responses for all model variables plus asset holdings and excess
% returns are stored in the cell array IPR. Each cell of IPR contains a
% cell array of impulse responses for each of the exogenous variables of
% the model.