%%%  Sectoral Reallocation Model : Phelan & Trejos (2000)  %%% 
%%  By Ryo Kato April, 2003 
%%  Phelan, C. and Trejos, A. (2000) in JME 45, 249-268.
%%   ``The Aggregate Effects of Sectoral Reallocations''

%%  This code is written by Ryo Kato (ryou.katou@boj.or.jp)
%%  You can use/rewrite this code without my permission.
%%  I will not be responsible for any damage/liability 
%%  induced by running this conde.

%%  This code is not for an impulse response of the system, but 
%%  to derive transition path from initial SS with old parameters 
%%  to the new SS with the new parameters. 
%%  Basically, what we do here on the code is to give the 
%%  old SS values of the state variables as the INITIAL values 
%%  to the new system with new parameters.
 clear all;

%%  Here is the shock you specify !!
  dd = -0.2 ;
 %% dd = NEWalpha - OLDAlpha = OLDalpha - NEWbeta 
 %%  that is, raising preference in good 1 and lowering that in good 2 
 %%  with unchanged preference for good 3


%% ------------------- [1] Parameter proc  ------------------------
%%  These are ``NEW'' (=post shock) parameters !! 

  alpha = 0.4;   %%  Preference param for Good1
  beta =  0.35;   %%  Preference param for Good2
  delta = 0.99;  %%  discount factor
  lamda = 2;     %% labor supply param
  myu = 2        %% labor supply elas > 1
  psi =0.3       %%  matching parameter
  gamma = 0.12  %% job quit rate 
  rho = -2;      %% elas of substitution = 1/(1-rho)  
  pg = 1 -alpha -beta;  %% Preference param for Good3
  
  %% Technology
  z1 = 0.8;
  z2 = 1.2;
  z3 = 1;
  
%% ---------------------<< [2] Steady State proc >> -----------------------
% SS shares of employment 
  G2 = ((alpha/beta)^(1/(rho-1)))*(z2/z1)^(rho/(rho-1))
  G3 = ((alpha/pg)^(1/(rho-1)))*(z3/z1)^(rho/(rho-1))

  G1 = 1 ;
  G = 1 + G2 + G3;
% SS hours and msearch efforts 
  hstar=(((alpha*((z1*G1)^rho)+beta*((z2*G2)^rho)+pg*(z3*G3)^rho)^(1/rho))/(lamda*myu*G))^(1/(myu-1));
  mstar = delta*(1-psi)*gamma*hstar*(myu-1)/myu;
% SS unemployment rate  
  ustar = ((mstar^psi)/gamma)^(1/(psi-1));
 % SS employment 
  N1star = (1-ustar)/G;
  N2star = G2*N1star;
  N3star = G3*N1star;
  Nstars = [N1star N2star N3star]';
  Nstar = sum(Nstars);
  SSvalues = [[ hstar  ; mstar ; ustar ] Nstars]
  
  
%% --------------------------[2] MODEL proc-----------------------------%%  
   % Define endogenous vars  ('a' denotes t+1 values)
    syms  h1  h2 h3  m1  m2 m3 N1 N2 N3 ha1 ha2  ha3  ma1 ma2 ma3 Na1 Na2 Na3; 
    syms  x1 x2 x3 xa1 xa2 xa3 ;
  %%  Definitions
  %% unemployment & Utility  
     ut = 1 - N1 - N2 -N3;
     ua = 1 - Na1 - Na2 -Na3;
     U = (alpha*x1^rho  + beta*x2^rho  + (1-alpha-beta)*x3^rho)^(1/rho) ;
     Ua= (alpha*xa1^rho + beta*xa2^rho + (1-alpha-beta)*xa3^rho)^(1/rho) ;
   %%  marginal utilities  
       U1 = diff(U,x1);
       U2 = diff(U,x2); 
       U3 = diff(U,x3);
       Ua1= diff(Ua,xa1);
       Ua2= diff(Ua,xa2);
       Ua3= diff(Ua,xa3);
   %%  Define production functions    
     x1s  = z1*N1*(h1 - m1);
     x2s  = z2*N2*(h2 - m2);
     x3s  = z3*N3*(h3 - m3);  
     xa1s = Na1*(ha1 - ma1);
     xa2s = Na2*(ha2 - ma2);
     xa3s = Na3*(ha3 - ma3);  
   %% Eliminating xi's from Ui's  
     U1q  = subs(U1,x1,x1s);
     U1p  = subs(U1q,x2,x2s);
     U1   = subs(U1p,x3,x3s);
     U2q  = subs(U2,x1,x1s);
     U2p  = subs(U2q,x2,x2s);
     U2   = subs(U2p,x3,x3s);
     U3q  = subs(U3,x1,x1s);
     U3p  = subs(U3q,x2,x2s);
     U3   = subs(U3p,x3,x3s);
     
     Ua1q  = subs(Ua1,xa1,xa1s);
     Ua1p  = subs(Ua1q,xa2,xa2s);
     Ua1   = subs(Ua1p,xa3,xa3s);
     Ua2q  = subs(Ua2,xa1,xa1s);
     Ua2p  = subs(Ua2q,xa2,xa2s);
     Ua2   = subs(Ua2p,xa3,xa3s);
     Ua3q  = subs(Ua3,xa1,xa1s);
     Ua3p  = subs(Ua3q,xa2,xa2s);
     Ua3   = subs(Ua3p,xa3,xa3s);
     
  %%% Optimal Conditions  & state transition %%%    
  %% Eliminatting marginal value functions
      V1  = (1/(delta*psi))*U1*(m1/ut)^(1-psi);
      V2  = (1/(delta*psi))*U2*(m2/ut)^(1-psi);
      V3  = (1/(delta*psi))*U3*(m3/ut)^(1-psi);
      Va1 = (1/(delta*psi))*Ua1*(ma1/ua)^(1-psi);
      Va2 = (1/(delta*psi))*Ua2*(ma2/ua)^(1-psi);
      Va3 = (1/(delta*psi))*Ua3*(ma3/ua)^(1-psi);
 %%%%%%%%%%%%%%%%%%%%%%%%%%% <Model Equations> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     
 %% Efficiency conditions & Envelope theorem & State transitions     
  oh1  =  z1*U1 - lamda*myu*h1^(myu-1);    
  oh2  =  z2*U2 - lamda*myu*h2^(myu-1);
  oh3  =  z3*U3 - lamda*myu*h3^(myu-1);   
  Value1   =  z1*U1*(h1-m1) - lamda*h1^myu ...
         + delta*Va1*((1-gamma) + (m1^psi)*(ut^(1-psi)))...
         - delta*(1-psi)*(ut^(-psi))*( Va1*m1^psi +Va2*m2^psi + Va3*m3^psi) - V1; 
  Value2   =  z2*U2*(h2-m2) - lamda*h2^myu ...
         + delta*Va2*((1-gamma) + (m2^psi)*(ut^(1-psi)))...
         - delta*(1-psi)*(ut^(-psi))*( Va1*m1^psi +Va2*m2^psi + Va3*m3^psi) -V2; 
  Value3   =  z3*U3*(h3-m3) - lamda*h3^myu ...
         + delta*Va3*((1-gamma) + (m3^psi)*(ut^(1-psi)))...
         - delta*(1-psi)*(ut^(-psi))*( Va1*m1^psi +Va2*m2^psi + Va3*m3^psi) - V3; 
 ST1 = N1*((1-gamma)+(m1^psi)*(ut^(1-psi))) - Na1;     
 ST2 = N2*((1-gamma)+(m2^psi)*(ut^(1-psi))) - Na2;     
 ST3 = N3*((1-gamma)+(m3^psi)*(ut^(1-psi))) - Na3;     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%  Compile equations      
      optcon = [ oh1 ; oh2 ; oh3 ; Value1 ; Value2  ; Value3  ; ST1 ; ST2 ; ST3 ];
      
%% ==========================  [3] Linearization proc ========================%%     
  % Differentiation 
     xx = [ha1 ha2  ha3  ma1 ma2 ma3 Na1 Na2 Na3 h1 h2 h3  m1  m2 m3 N1 N2 N3 ] ; % define vars
     jopt = jacobian(optcon,xx);
     
   % Evaluate each derivative
     ma1 = mstar;  ma2 = mstar; ma3 = mstar;
     m1  = mstar;  m2 = mstar;  m3  = mstar;
     h1  = hstar;  h2 = hstar;  h3  = hstar;
     ha1 = hstar;  ha2 = hstar; ha3 = hstar;
     Na1 = N1star;  Na2 = N2star;  Na3 = N3star;
     N1  = N1star;  N2 = N2star;    N3  = N3star;
     
 % Define Linear Coefficients  
     coef = eval(jopt);
   % In terms of % deviations from ss    
    vo = [ mstar  mstar  mstar  hstar  hstar hstar  N1star N2star N3star  ] ;
    TW = [ vo ; vo  ; vo ; vo  ; vo ; vo  ;  vo  ; vo ; vo] ;
 
    B =  [ -coef(:,1:9)  ].*TW ;
    C =  [  coef(:,10:18)  ].*TW ;
 % B[c(t+1)  l(t+1)  k(t+1)  z(t+1)] = C[c(t)  l(t)  k(t)  z(t)]
    A = inv(C)*B  ;%(Linearized reduced form )
   
%% =========== [4] Solution proc (Do NOT touch!!) ============== %%
%  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))<0.9999999999
         	UQ(jjw,:) = Q(jj,:);
         	jjw = jjw+1;
      	end
      end
    % Extract stable roots 
      VLL = [];    jjjw = 1;
      for jjj = 1:length(theta)
      	if abs(theta(jjj)) >1.0000000001
            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)
   VL = inv(diag(VLL))
   
   % 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)   ;
    
%% =================== [5]  SIMULATION proc ========================== %%    
%  TIME&INITIAL VALUES 
     t = 24;			% Time span   
     
 %%%%%%%%%%%  Compute the OLD SS for initial values %%%%%%%%%    
%%  These are ``OLD' (=pre-shock) parameters !! 
%%  ======== Here are the SHOCKS you specify!! ======== %%
  alphaX = alpha - dd ;   %%  Preference param for Good1 + X
  betaX  = beta  + dd ;   %%  Preference param for Good2 - X
  gammaX = gamma  %% job quit rate 
  rhoX   = rho;      %% elas of substitution = 1/(1-rho)  
  pgX    = 1 -alphaX -betaX;  %% Preference param for Good3
  
  %% -----------------<< OLD Steady State proc >> --------------  %%
  % SS shares of employment   
   G2X = ((alphaX/betaX)^(1/(rhoX-1)))*(z2/z1)^(rho/(rho-1));
   G3X = ((alphaX/pgX)^(1/(rhoX-1)))*(z3/z1)^(rho/(rho-1));
   G1X = 1 ;
   GX = 1 + G2X + G3X;
 % SS hours and msearch efforts 
   hstarX=(((alphaX*((z1*G1)^rhoX)+betaX*((z2*G2)^rhoX)+pgX*(z3*G3)^rhoX)^(1/rhoX))/(lamda*myu*G))^(1/(myu-1));
   mstarX = delta*(1-psi)*gammaX*hstarX*(myu-1)/myu
 % SS unemployment rate  
   ustarX = ((mstarX^psi)/gammaX)^(1/(psi-1))
  % SS employment 
   N1starX = (1-ustarX)/GX;
   N2starX = G2X*N1starX;
   N3starX = G3X*N1starX;
   NstarX = N1starX + N2starX + N3starX;
 %%  ===================================================== %%   
 %%  Converting initial SS in terms of deviations from the NEW SS.
  NS1 = (N1starX - N1star)/N1star;
  NS2 = (N2starX - N2star)/N2star;
  NS3 = (N3starX - N3star)/N3star;
   
 %% Initial values of state varialbe   
   S1 = [ NS1   NS2   NS3 ]' ;  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
% SIMULATION/Impulse responses 
      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')';
      
%%  Reverting levels
      
  Stars =  ones(length(SY),1)*Nstars'   ;   
  LSY   = SY.*Stars + Stars
  N = sum(LSY')'
  
  Stars2 = ones(length(X),6).*hstar ;
  HV = X.*Stars2 + Stars2;
  HX = HV(:,1:3).*LSY    %%  labor*hours   
  HLX = sum(HX')'
  AHs = HLX./sum(LSY')' ;
  AH = [hstarX ; AHs(1:t) ];
  
  HL = [hstarX*NstarX ; HLX(1:t)]
  
  Hi = [hstarX hstarX hstarX  ; HV(1:t,1:3)]
  
%% ===============================  
  
%  DRAWING FIGURES
   figure(1)
    subplot(3,1,1);
       plot( LSY );
       ylabel('Ni')
       xlabel('time');
    subplot(3, 1,2);
       plot(AH);
       ylabel('hi')
       xlabel('time');
    subplot(3, 1,3);
       plot(Hi);
       ylabel('hi')
       xlabel('time');  

 S2 =  [ (N1star-N1starX)/N1starX (N2star-N2starX)/N2starX  (N3star-N3starX)/N3starX ]







   
      
