%%  Staggered Wage model %% 
%%  By Ryo Kato and Takayuki Tsuruga 
%%  October 30, 2002 

%%  This code is written by Ryo Kato (kato.13@osu.edu)
%%  & Taka Tsuruga (tsuruga.1@osu.edu)
%%  You can use/rewrite this code without our permission.
%%  I will not be responsible for any damage/liability 
%%  induced by running this conde.

clear all;

%% -----------------------------------------------------------
% [1]  Parameters
  sigma = 1;   % CRRA consumption 
  phi   = 2;   % CRRA labor
  e     = 10  ;  %  Goods markt elasticity 
  ew    = 2   ;  %  Labor makt elasticity. 
  c     = log(e/(e-1));   
  cw    = log(ew/(ew-1)); 
  beta    = 0.99;    % discount rate
  rho     = 0.8;     % AR(1) in technology shock
  rhom    = 0.8;     % AR(1) in money level
  delta   = 0.25 ;   % frequency of p change
  deltaw  = 0.2  ;   % frequency of w change
  gamma   = (1-beta*(1-delta))*delta/(1-delta) ;  
  gammaw  = (1-beta*(1-deltaw))*deltaw/(1-deltaw) ;

%% --------------------- << SS value proc >> -----------------------
     mstar   = 0 ;
     vstar   = 0 ;
     pstar   = -(c+cw)/(phi+sigma) ;
     ystar   = -pstar ;
     lstar   = ystar ;
     wstar   = pstar - c;
     mcstar  = -c;
     mrsstar = -c-cw;
     
    YWoMCMRS = [ ystar  wstar ; mcstar  mrsstar ]  % show SS values
     
%% -------------- << Linearization proc  >> ---------------------------
%% ---------------------------- MODEL -----------------------------%%  
   % Define endogenous vars  ('a' denotes t+1 values)
   syms  xt  xwt  lt  pt  wt   mt  vt ... 
         xa  xwa  la  pa  wa   ma  va     ; 
   % MC & MRS  
       mca  = wa - pa - va ;      % marginal cost(t+1)
       mrsa = phi*la + sigma*(va + la) ;  % marginal rate of s(t+1)
       
   % Optimal Conditions  & state transition
       pai  = (pa - pt) - gamma*(mca + c)    - beta*(xa-xt) ;   % NKPC    
       dw   = (wa - wt) - gammaw*(mrsa + cw - wa + pa ) - beta*(xwa-xwt) ; % Wage Phill curve
       dump  = pa - xt  ;  % p(t+1) dummy for companion form
       dumw  = wa - xwt ;  % w(t+1)  dummy for companion form
       labor = lt + vt - mt + pt ;  % M-P = Y  where Y = vL
       money = ma - rhom*mt  ; % Monetary policy
       tech  = va - rho*vt   ; % technology shock       
       
       optcon  = [pai ; dw  ; dump  ; dumw ; labor ; money ; tech ];    
       
%% ----------------------------------------------------------------%%     
  % Differentiation 
     xx = [xa  xwa  la  pa  wa ma  va  xt xwt lt pt wt  mt vt  ] ; % define vars
     jopt = jacobian(optcon,xx);

   % Evaluate each derivative
     xa  = pstar;     xt  = pstar;
     xwa = wstar;     xwt = wstar;
     la  = lstar;     lt  = lstar;
     pa  = pstar;     pt  = pstar;
     wa = wstar ;     wt = wstar ;
     ma = mstar ;     mt = mstar ;
     vt = vstar;   va = vstar;
     
   % Define Linear Coefficients  
     coef = eval(jopt);
     
   % In terms of % deviations from ss    
    vo = [ 1  1  1  1  1  1  1 ] ;
    TW = [ vo ; vo  ; vo ; vo  ; vo  ;vo ; vo ] ;
 
    B =  [ -coef(:,1:7)   ].*TW ;
    C =  [  coef(:,8:14)  ].*TW ;
  % B[p(t+2)  w(t+2)  l(t+1)  p(t+1)  w(t+1)  m(t+1)  v(t+1)] 
  % = C[p(t+1)  w(t+1)  l(t)  p(t)  w(t)  m(t)  v(t)]
    A = inv(C)*B;  %(Linearized reduced form )
   
%% ================== <<Solution proc>> ======== %%
% [2] EIGEN DECOMPOSITION
   [W V] = eig(A);
   Q = inv(W);
   W*V*Q; 
   theta = diag(V)
  % Extract stable vectors
      SQ = [];    jw = 1;
      for j = 1:length(theta)
      	if abs(theta(j)) >1.000000001
         	SQ(jw,:) = Q(j,:);
         	jw = jw+1;
      	end
      end
   % Extract unstable vectors
     UQ = [];   jjw = 1;
      for jj = 1:length(theta)
      	if abs(theta(jj))<1.000001
         	UQ(jjw,:) = Q(jj,:);
         	jjw = jjw+1;
      	end
      end
    % Extract stable roots (lamda>1) note GYAKU!
      VLL = [];    jjjw = 1;
      for jjj = 1:length(theta)
      	if abs(theta(jjj)) >1.000000001
            VLL(jjjw,:) = theta(jjj,:); 
            jjjw = jjjw+1;
      	end
      end
      
   % Show Eigen Vectors on U-S Roots  
     UQ;  % n x n+k
     SQ;  % k x n+k 
     
  % [3] ELIMINATING UNSTABLE VECTORS
      k = min(size(SQ));   % # of predetermined vars
      n = min(size(UQ));   % # of jump vars
      nk = [n k]
   % Stable V (eig mat)
      VLLL = [];   hh = 1;
   for h = 1:k
      VLLL(:,hh) = VLL(h)*ones(k,1);
      hh = hh +1 ;
   end
   VL = inv(eye(k).*VLLL) ; % modified inverse!!
   
   % Elements in Q
    PA = UQ(1:n,1:n);    PB = UQ(1:n,n+1:n+k);
    PC = SQ(1:k,1:n);    PD = SQ(1:k,n+1:n+k);
    P = -inv(PA)*PB ; % X(t) = P*S(t)
    PE = PC*P+PD ;
   
   % SOLUTION
    PX = inv(PE)*VL*PE;
    AA = real(PX)   

    
%% ----------<< SIMULATION PROGRAM  >>------------ %%    
% [4] TIME&INITIAL VALUES 

     t = 24;			% Time span you want to simulate   
    
   % Initial Values 
   % state var = [p(t)  w(t)  m(t)  v(t)]
   S1 = [ 0   0  0.1   0]' ;

    
% [5] SIMULATION 
      Ss = S1;
      S = zeros(t,k) ;  
        for i = 1:t
          q = AA*Ss  ;
   	      S(i,:) = q';
          Ss = S(i,:)';
     	end
      SY = [S1' ;S]  ; 
      X = (real(P)*SY')';
      
      infl  = X(:,1)- SY(:,1) ;  % define inflation
      dwage = X(:,2)- SY(:,2) ;  % define wage growth
      
% [6] DRAWING FIGURES
   figure(1)
     subplot(2,1,1);
       plot( [X(:,3)  SY(:,4)+X(:,3) ]);
       ylabel('Lt(bl)  Yt(gr)')
       xlabel('time');
     subplot(2, 1,2);
       plot( [infl  dwage]  );
       ylabel('dp(blue)  dw(grn)')
       xlabel('time');
     



