% oligopoly2
% Sargent's amendment of oligopoly.m, Stijn's original program.
% See comments signed by tjs below and the new section
% The approximating model is computed.  See comments below.
% See oligop2.ps for description
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This program computes the oligopoly example
% Robustness in forward looking models
% Stijn Van Nieuwerburgh
% First version: January 27, 2001 (stijn)
% Current version: January 31, 2001 (sargent)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% There is a representative competitive firm with cost function
%  d q_t + .5 h q_t^2 + .5 c(q_{t+1} - q_t)^2 
% A large firm has cost function
%  eQ_t + .5 g Q_t^2 + .5 c(Q_{t+1} - Q_t)^2
% Demand is 
% p(t)= ao-a1*(Q(t)+q(t))+v(t)
% v(t+1)= rho*v(t)+Ce*eps(t+1); eps is iid(0,1)
%
% Equilibrium concept: Stackelberg, with mu_{xt} being the implementation
% multiplier. 
% the state vector in the final representation is
% y_t' = [1 v_t Q_t q_t \mu_{xt}]' 

     a0= 100;    % intercept of demand function
     a1=1;      % slope of demand function
     rho=0.8;   % autocorrelation factor
     Ce=.2;     % variance of shocks
     c=1;      % adjustment cost factor
     b=0.95;    % discount factor
     d=20;      % unit cost for competitive firm
     e=20;      % unit cost for large firm
     g=.2;     %  quadratic term in cost function for large firm
     h=.2      % quadratic term in cost function for competitive firm
     
     th=10000000;     % robustness parameter
    th=10;
    %th=10000000;
% Define the following matrices
     AA1=[1 0 0 0 0 ;0 1 0 0 0 ;0 0 1 0 0;0 0 0 1 0;(a0-d) 1 -a1 -a1-h c];
     AA2=[1 0 0 0 0;0 rho 0 0 0;0 0 1 0 0;0 0 0 1 1 ;0 0 0 0 (c/b)];
        % tjs changed last entry from b/c to c/b; see the notes
     BB=[0;0;1;0;0];
     CC=[0;Ce;0;0;0];
     A=inv(AA1)*AA2;
     Bb= inv(AA1)*BB;
     C = inv(AA1)*CC;
     B=inv(AA1)*[BB CC];
     
% make R and Q positive definite (in absence of be*th 
     R=[c/2 0;0 -b*th] ;  
     Q= -[0 0 ((a0-e)/2) 0 0;0 0 .5 0 0;((a0-e)/2) .5 -a1-g/2 -(a1/2) 0; 0 0 -(a1/2) 0 0;0 0 0 0 0];
               % minus sign added by tjs
               
     
% Define the number of 'jump variables'.
     nx=1;
     n= max(size(Q));
     nz=n-nx;

% Solve the optimal linear regulator from the above matrices
     [F P]= olrp(b,A,B,Q,R);

% Manipulate matrices

P21=P(nz+1:n,1:nz);
P22=P(nz+1:n,nz+1:n);
F1 =F(1,:) ;
F2 =F(2,:);
A1 =A(1:nz,:);
A2 =A(nz+1:n,:);
B1 =Bb(1:nz,:);
B2 =Bb(nz+1:n,:);
C1 =C(1:nz,:);
C2 =C(nz+1:n,:);

% Form the matrices M and MM

MW=[eye(nz) zeros(nz,n-nz); P21 P22]*[A1-B1*F1-C1*F2 ;A2-B2*F1-C2*F2]*[eye(nz) zeros(nz,n-nz);-inv(P22)*P21 , inv(P22)];
M=[eye(nz) zeros(nz,n-nz); P21 P22]*[A1-B1*F1 ;A2-B2*F1]*[eye(nz) zeros(nz,n-nz);-inv(P22)*P21 , inv(P22)];
%  MW is the worst-case law of motion, corrected from Stijn's program
%  M is the approximating law of motion.  tjs 1/30/01
%
%  Below, we'll use M to simulate the approximating model under the
%  robust rule.  We'll use Mw to verify an identity satisfied by
%  decisions $i_t$ for the representative firm

MM= [-inv(P22)*P21 , inv(P22)];
MM11=M(1:nz,1:nz);
MM12=M(1:nz,nz+1:n);
MM21=M(nz+1:n,1:nz);
MM22=M(nz+1:n,nz+1:n);

% Compute the ARMA representations for controls and perturbations
FF= -F*[eye(nz) zeros(nz,n-nz);-inv(P22)*P21 , inv(P22)];
FF11=FF(1,1:nz);
FF12=FF(1,nz+1:n);
FF21=FF(2,1:nz);
FF22=FF(2,nz+1:n);

  %%%%  The following  formulas altered by tjs (use pinv)
  %  We can use pinv function (generalized inverse)
  % I've %'d out Stijn's formulas
  
%H1=inv(FF12'*FF12)*FF12';
%H2= FF12*H1;
%H3=inv(FF22'*FF22)*FF22';
%H4= FF22*H3;

%AR11=inv(H2)*FF12*MM22*H1;
%MA11=FF11;
%MA12=inv(H2)*FF12*(MM21-MM22*H1*FF11);

%AR21=inv(H4)*FF22*MM22*H3;
%MA21=FF21;
%MA22=inv(H4)*FF22*(MM21-MM22*H3*FF21);
 
AR11=FF12*MM22*pinv(FF12);
MA11=FF11;
MA12=FF12*(MM21-MM22*pinv(FF12)*FF11);



AR21=FF22*MM22*pinv(FF12);         %  formulas altered by tjs
MA21=FF21;
MA22=FF22*(MM21-MM22*pinv(FF12)*FF11);




disp('The feedback representation for u is:');
disp('u(t)=');disp(FF11);disp('* z(t) + ');
disp(FF12);disp('*mu_x(t)');

disp('The feedback representation for w is:');
disp('u(t)=');disp(FF21);disp('* z(t) + ');
disp(FF22);disp('*mu_x(t)');




disp('The ARMA representation for u is:');
disp('u(t)=');disp(AR11);disp('* u(t-1) + ');disp(MA11);disp('* z(t) + ');disp(MA12); disp('* z(t-1) .');




disp('The ARMA representation for w is:');
disp('w(t)=');disp(AR21);disp('* u(t-1) + ');disp(MA21);disp('* z(t) + ');disp(MA22); disp('* z(t-1) .');

 % above corrected by tjs  

%%%%%% new part of program 
 
% Following material added by tjs, 01/30/01
%% The new stuff consists of two parts.  

% Part I.  Use sys to form state space system.  Then get impulse
% response and simulation of system.

% first form the volatility matrix for z_t, mu_{xt}  by premultiplying

Ct=[eye(nz) zeros(nz,n-nz); P21 P22]*C;
%
%  (this transformation follows from the transformation 
%  from [z_t, x_t] to [z_t, \mu_{xt}]  .  See chapter 10.  
%
%  get state space system for p, q, Q, mu_x, w, v
%
G=[a0 1 -a1 -a1 0 ;    % picks off p
   0 0 0 1 0 ;         % selects q
   0 0 1 0 0];          % selects Q
G=[G; [0 0 0 0 1]];     % selects \mu_x 
%   G=[G; FF11 FF12];   % last row picks off u=(Q_{t+1} - Q_t)
G=[G; FF21 FF22];       % picks off w_{t+1}, worst case shock
G=[G;0 1 0 0 0];        % picks off v_t, the demand shock
% now form the state space system M, Ct, G, DD=zeros(6,1)


DD=zeros(6,1);
sys1=ss(M,Ct,G,DD,1) ;   % system for approximating model
sys2=ss(MW,Ct,G,DD,1);   % worst case model
 
%  now compute impulse response functions for p, q, Q, u with respect
% to the scalar demand shock epsilon.  

nn=12;   %  mm is length of impulse response
y1=impulse(sys1,nn);  y1=y1(2:nn,:);
y2=impulse(sys2,nn);  y2=y2(2:nn,:);
figure(1)

tt=[0:nn-2];  % for cosmetics of figure, to label x-axis properly

figure(1)   %  impulse response
subplot(4,1,1)
stairs(tt,y1(:,1)),ylabel('p')
subplot(4,1,2)
stairs(tt,y1(:,2)),ylabel('q')
subplot(4,1,3)
stairs(tt,y1(:,3)),ylabel('Q')
subplot(4,1,4)
stairs(tt,y1(:,4)),ylabel('\mu_x')

print -depsc oli1

figure(6)  % rest of impulse response
subplot(3,1,1)
stairs(tt,y1(:,2)+y1(:,3)),ylabel('q +Q')
subplot(3,1,2)
stairs(tt,y1(:,5)),ylabel('w')
subplot(3,1,3)
stairs(tt,y1(:,6)),ylabel('v')
 
print -depsc oli6

%  Now compute simulation, starting from appropriate initial state.



TT=50+1;
u=randn(TT,1);

tp=[0:TT-1];
%  construct steady state

xs=null(eye(5)-M);
xs=xs./xs(1);
%    
%   xs is steady state value of  [1 v Q q \mu_x];  

disp('theta = ' ),th
disp('steady state of price, q, Q')

ps=[a0 1 -a1 -a1  0]*xs
qs=[0 0 0 1 0]*xs
Qs=[0 0 1 0 0]*xs

disp('steady state xs'),xs

xinitial=xs; xinitial(5)=0;   % start from zero shadow price; 
                              %  all other state variables from steady state
                              
%  now simulate and prepare plots

Y1=lsim(sys1,u,[],xinitial);

figure(2)
subplot(3,1,1)
plot(tp,Y1(:,2)+Y1(:,3),'-g'),ylabel('q+Q')
subplot(3,1,2),
plot(tp,Y1(:,2),'-b'),ylabel('q')
subplot(3,1,3)
plot(tp,Y1(:,3),'--r'),ylabel('Q')
print -depsc oli2

figure(3)
subplot(3,1,1)
plot(tp,Y1(:,4),'-g'),ylabel('\mu_x')
subplot(3,1,2)
plot(tp,Y1(:,3),'-b'),ylabel('Q')
subplot(3,1,3)
plot(tp,Y1(:,1),'--r'),ylabel('p')
print -depsc oli3


figure(4)
subplot(2,1,1)
plot(tp,Y1(:,5)+Y1(:,6),'--r');ylabel('w +v')
hold on
plot(tp,Y1(:,6),'-b');ylabel('v')
hold off
subplot(2,1,2)
plot(tp,Y1(:,5),'--r');ylabel('w')

print -depsc oli4
figure(5)
subplot(4,1,1)
plot(tp,Y1(:,6),'-g');ylabel('v');%title('path of v')
subplot(4,1,2)
plot(tp,Y1(:,3),'-b'),ylabel('Q')
subplot(4,1,3)
plot(tp,Y1(:,2),'-b'),ylabel('q')
subplot(4,1,4)
plot(tp,Y1(:,1),'--r'),ylabel('p')

print -depsc oli5


%%  Part II.  Perform check on decision rule, mentioned in
%  chapter 10.


%  first get the decision rule for i_t
% i_t = - P22^{-1} P21 z_t + P_{22}^{-1} \mu_{xt}

if1= -inv(P22)*P21
if2=inv(P22)

iff=[if1 if2]

%% now get roots of characteristic polynomial

p=[1 -(b^(-1)+1+c*h)  b^(-1)]
rr=roots(p)
lam=rr(2)
  %  note: rr(2) = lam; 1/rr(1) = b*lam
  
%  form vector to pick off p-d from [z mu_x];

ep=[a0-d 1 -a1 -a1 0];

 % now form decision rule by using forecasts

ig1=[ 0 0 0 (lam-1) 0];

ig2=c^(-1)*ep*b*lam*MW*inv(eye(5) - b*lam*MW);

igg=ig1+ig2

iff
  
disp('The following should be zero'), 
iff-igg
   %  igg should agree with iff
return
