%-----------------------------------------------------------------------------
%
% Program: test1.m
%
%  Solves Svensson's 1997 model with robustness.
%  The model is given by the two equations
%
%  x(t+1)  = b1x(t)-b2(i(t)-pi(t))+u(t+1)
%  pi(t+1) = pi(t)+a1x(t)+e(t+1)
%
%  where x is the output gap, i the nominal interest rate, pi the inflaiton rate,
%  u an AD shock and e a cost push shock.
%
%  The central banker sets i. The loss function for the central banker is
%  quadratic in x , pi, i and i(t)-i(t-1).
%
%  Paolo Giordani and Paul Söderlind, February 2002
%-----------------------------------------------------------------------------

theta = 10;        %degree of robustness. 5-100
AD = 1;            %1 for AD shock, else cost-push shock

Tbig  = 10;        %Periods to calculate impulse response for

%coefficients of the law of motion

a1 = 0.4;          %Ball (1997) uses a1=0.4 b1=0.8 b2=1 coeffinfl=1
b1 = 0.8;
b2 = 1;
coeffinfl = 1;    %the coefficient on lagged inflation in the Phillips curve

%coefficients for matrix C

stdinfl = 0.5;     %std of the cost-push shock
stdy = 0.7;        %of AD shock
stdi = 0.7;        %of the ad hoc monetary policy shock

%coefificents in loss function

bet = 0.98;          %discount factor
outputweight = 0.5;  %weight of squared output gap in the loss function
drweight = 0.2;      %coefficient on (i(t)-i(t-1))^2 in loss function
rweight = 0;         %coefficient on i(t)^2 in loss function

cutoff = 1.01;
%----------------------------------------------------------------------------

%writing the state-space form
[A,B,Q,U,R,uy,upi,C] = Svensson(a1,b1,b2,coeffinfl,stdinfl,stdy,outputweight,rweight,drweight);

if AD == 1;
  shock=C*uy;
else;
  shock=C*upi;
end;
%----------------------------------------------------------------------------

n1=3; n2=0;  %no forward-looking variable

[M,N,Ma,Na,Fv,J0,J0a] = ComAlgSR( A,B,Q,R,U,bet,cutoff,C,theta );

x1w       = Var1SimPs(M,shock,Tbig);
interestw = x1w*N';
interestw = interestw(:,1);
Simupi_Dw = [x1w(:,1:2),interestw];   %IRF of pi(t),y(t),i(t) to shock, worst case
%----------------------------------------------------------------------------

x1a       = Var1SimPs(Ma,shock,Tbig);         % approximating model
interest  = x1a*N';
interest  = interest(:,1);
Simupi_Da = [x1a(:,1:2),interest];             % IRF of pi(t),y(t),i(t) to price shock


%Plotting graph
xx = (1:Tbig)';
if AD == 1;
  aux = 'AD shock';
else;
  aux = 'CP shock';
end;


if exist('OCTAVE_VERSION');       %if Octave
  gnuplot_has_multiplot = 1;
  automatic_replot = 0;
  multiplot(2,2);
  gset('nokey');
  title('');
    title(['Output gap to ' aux]);
    mplot(xx,[Simupi_Da(:,1),Simupi_Dw(:,1)]);
    title(['Inflation to ' aux]);
    mplot(xx,[Simupi_Da(:,2),Simupi_Dw(:,2)]);
    title(['Short rate to ' aux]);
    mplot(xx,[Simupi_Da(:,3),Simupi_Dw(:,3)]);
    multiplot(0,0);
    gset('key');
    automatic_replot = 1;
else;                                  %if MatLab
  figure;
    subplot(2,2,1);
    plot(xx,[Simupi_Da(:,1),Simupi_Dw(:,1)]);
    title(['Output gap to ' aux]);
    subplot(2,2,2);
    plot(xx,[Simupi_Da(:,2),Simupi_Dw(:,2)]);
    title(['Inflation to ' aux]);
    subplot(2,2,3);
    plot(xx,[Simupi_Da(:,3),Simupi_Dw(:,3)]);
    title(['Short rate to ' aux]);
end;
%----------------------------------------------------------------------------