%--------------------------------------------------------------------------
%  sosdebtbeg.m
%
%  Program uses functions in ComItAlg.m, DiscAlg, and SimpRulT.m.
%  These functions call on qzswitch.m (by Simms) and reorder.m (by Klein).
%  The program also uses Var1SimPs to calculate impulse response functions.
%
%  The program uses the model in Goyal, A., 2009, "Sustainable debt and deficits in emerging markets". 
%  It adds a debt equation to forward looking dynamic aggregate demand and supply equations, 
%  calibrated for the Indian economy.
%  
%  It was modified from a program originally written by Paul Soderlind (nepso@hhs.se, 20 June 2000) 
%  for a simple version of the model in Fuhrer's 1997 JMCB paper.  
%
%  Ashima Goyal (ashima@igidr.ac.in, 13 July 2009) 
%----------------------------------------------------------------------------

diary('sosdebt.out');            %printing: to screen AND the file
%-----------------------------------------------------------------------

                      
al   = 0.5;              
eta  = 0.4;
phi  = 0.25;
ke   = 0.1;
lam  = 0.24;
bet  = 0.99;
rho  = 0.01;
gr   = 0.006;
tb   = 0.005;
bb   = 0.8;

sigr = 1;
%sig = sigr/eta
w = sigr+(1-al)*(sigr-1);
%siget = (sig*eta*(1-al) +al*sigr)/(eta*(1-al)+w*al)
sigd = sigr/(eta*(1-al)+w*al);
%sigal = sigr/((1-al)+w*al)
%mcs = (sig*eta*(1-al)+al*sigr)/sigr
%d = siget/(phi + siget)
%sum = d*((1-eta)*(sig-siget))

vard = 0.19^2;                   %dd shocks
vars  = 0.2^2;                    %cost shocks
vare = 0.2^2;
varg = 0.1^2;
varb = 0.2^2;

qy    = 0.7;                      %Loss function
qpi   = 2;                      %at=1; cp=0.08; y*=1
qf    = 1;


[A,B,K,Q,U,R] = sosdebt(al,eta,phi,ke,lam,bet,rho,gr,tb,bb,vard,vars,vare,varg,varb,qy,qpi,qf);

n1 = 8;
n2 = 4;
k  = 1;

Tbig   = 12;                       %periods to simulate in VAR


%Shock0 = [0;sqrt(vard);zeros(6,1)];   %Impulse response wrt ep

 %Shock0 = [0;0;0;sqrt(vare);zeros(4,1)]; %Impulse response wrt e shock

 %Shock0 = [0;0;sqrt(vars);0;zeros(4,1)];   %Impulse response wrt cost shock
 
 %Shock0 = [0;0;0;0;zeros(2,1);sqrt(varg);0];   %Impulse response wrt g exp shock

 Shock0 = [1;0;0;0;zeros(4,1)];            %Impulse response wrt rn shock
 
 %Shock0 = [zeros(7,1);sqrt(varb)];   %Impulse response wrt debt bt shock
%-----------------------------------------------------------------------
 
                             
                                  %DISCRETION

disp('Please wait, discretionary case takes some time to solve');
if exist('OCTAVE_VERSION');
  page_output_immediately = 1;   %seems to make Octave work better
end;
[M_Disc,C_Disc,V,F] = DiscAlg( A,B,Q,R,U,bet,n1,n2,eye(n1),zeros(n2,n1),...
                               [1E-1;1E-5],0,1,0,1,0,1e+4 );
F
x1 = Var1SimPs( M_Disc,Shock0,Tbig );  %VAR of x1(t)
x2 = x1*C_Disc';                        %x2(t)
x  = [x1, x2];                          %x(t) = [x1(t),x2(t)]
uu = - x1*F';                            %u(t)

%picypiie_Disc = [ x(:,5),x(:,7), x(:,8), x(:,9), uu ];

picypiie_Disc = [ ((1-al)*(x(:,10))+ al*(x(:,11))-al*(x(:,4))),x(:,9), x(:,10), x(:,12), uu ]; %fl cpi

tid = linspace(1,Tbig,Tbig)';
if exist('OCTAVE_VERSION');
  title('Discretion: impulse response to a one std of price shock');
  plot( tid,picypiie_Disc(:,1),'-1;Consumer Inflation;',...
        tid,picypiie_Disc(:,2),'-2;Output;',...
        tid,picypiie_Disc(:,2),'-3;Domestic Inflation;',...
        tid,picypiie_Disc(:,2),'-4;Exchange rate;',...
        tid,picypiie_Disc(:,3),'-5;Interest rate;');
else;
  figure(3);
  %plot( tid,picypiie_Disc );
  hold on
  subplot(3,1,1);
  %plot( tid,x(:,5),'-*k');
  plot( tid,((1-al)*(x(:,10))+ al*(x(:,11))-al*(x(:,4))),'-*k'); %does not report exchnage rate
  title('Figure: Benchmark DIT: cost shock');  
  hold on 
  plot( tid,x(:,9),'-+r');
  hold on 
  plot( tid,x(:,10),'--dm');
  hold on 
  plot( tid,x(:,12),':hb'); hold on 
  plot( tid,uu,'-^c');
  %title('Full "DIT": impulse response to a calibrated cost shock');
  legend('Consumer inflation','Output','Domestic inflation','G debt','Interest rate',-1);
 end;


disp(' ');disp('debt benchmark, +.o1 rn shock (-ve shock-for nat rate fall) with added -.01 term, rm -bb*Aj-bb*gr in debt eq and tb =.005(pi-c,y,pi,b,i)'); %- .01Aa impulse response to an rn shock:
disp( [linspace(1,Tbig,Tbig)', picypiie_Disc] );
Y = std(picypiie_Disc,1)
disp Y

%-----------------------------------------------------------------------

diary off;
