%% Carlstrom-Fuerst Model (Infinite horizon Bernanke-Gertler model)
%%  AER (1997)
%%  Ryo Kato March 18, 2002 

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

% Uniform distribution [0 2]

clear all;
close all;

%% -----------------------------------------------------------
% [1]  Parameters
  alpha  = 0.3;
  alphah = 0.01;
  beta   = 0.99;
  phi    = 0.9; 
  sigma  = 1.5;
  myu    = 0.2;
  delta  = 0.025;
  gamma  = 0.95;
  Lstar  = 1.5 ;
  
  ehta = 0.33 ; 

%% ----------- <<Optimal Contract>> ------------------- %%
   syms w x qt qstarsym  

   pdf = 1/2  ;  %pdf uniform [0 2]
   cdf = int(pdf,0,x); 
   fxs = int(pdf*w,w,x,2)-(1-cdf)*x  ;
   gxs = int(pdf*w,w,0,x)-cdf*myu+(1-cdf)*x;

   contract = qt - 1/(1-cdf*myu + (x/2)*myu*(fxs/diff(fxs))); % qt - h(x) = 0
   xs = solve(contract,x);
   x = xs(2)  % x = h(qt)

 fx  = eval(fxs)  % E's share f(x(qt))
 gx  = eval(gxs)  % L's share g(x(qt))
 cdf = eval(cdf)  % PHI(x(qt)) 
 
 rho = (qt*fx)/(1-qt*gx)-1
 rhoq = subs(rho,qt,qstarsym)
 
%% =========================== << SS value proc >> =========================== %%
%-------------------------------------------------------% 
%qstar from optimal contract
     ssq = 1-gamma*(1+rhoq); % q 
     qqs = solve(ssq,qstarsym);
     qstar = double(qqs(1));
%------------------------------------------------------ %  
     xstar   = subs(x,qt,qstar);
     cdfstar = subs(cdf,qt,qstar) ;
     fxstar  = subs(fx,qt,qstar);
     gxstar  = subs(gx,qt,qstar);
      
  rstar   = qstar*((1-beta/beta)+delta);
  kstar   =  double(((rstar/alpha)*Lstar^(alpha+alphah-1))^(1/(alpha-1)) );
  wstar   =  double( (1-alpha-alphah)*(kstar^alpha)*(Lstar^(-alpha-alphah))  );
  Ystar   =  double( (kstar^alpha)*(Lstar^(1-alpha-alphah)));
  rhostar = ((qstar*fxstar)/(1-qstar*gxstar))-1;
  istar   = delta*kstar/(1-myu*cdfstar);
  nstar   = (1-qstar*gxstar)*istar;
  Zstar   = (nstar-alphah*Ystar)/(qstar*(1-delta)+rstar);
  cestar  = (1+rhostar)*nstar-qstar*Zstar;
  cstar   =  Ystar-(istar+cestar)*ehta;
  lamda  = wstar*cstar^(-sigma) ;
  
  QFGvRKR2vINZvCC = [qstar  fxstar  gxstar  ; rstar   kstar  rhostar ; 
                     istar  nstar   Zstar   ; cestar  cstar  xstar   ;
                     Ystar  lamda   Lstar]'
                     
%% ============================= << MODEL >> ================================ %%  
 
 syms  ct Kt vt qt ca Ka va qa Zt Za  
                                              
  %% Preparation
    fxt  = fx;
    gxt  = gx;
    fxa  = subs(fxt,qt,qa) ;
    gxa  = subs(gxt,qt,qa) ;
    cdft = cdf;
    
    %%  Market Clear 
  Lt = ((1-alpha-alphah)*vt*(Kt^alpha)/(lamda*(ct^sigma)))^(1/(alpha+alphah)) ;
  La = ((1-alpha-alphah)*va*(Ka^alpha)/(lamda*(ca^sigma)))^(1/(alpha+alphah)) ;
  rt = alpha*vt*(Kt^(alpha-1))*(Lt^(1-alpha-alphah));
  ra = alpha*va*(Ka^(alpha-1))*(La^(1-alpha-alphah));
  Yt = vt*(Kt^alpha)*(Lt^(1-alpha-alphah));
  nt = alphah*Yt + Zt*(qt*(1-delta)+rt);  
  it = (1/(1-qt*gxt))*nt;
  cet =  ((Yt - ct)/ehta) - it;
  rhot = (qt*fxt/(1-qt*gxt))-1;
  rhoa = (qa*fxa/(1-qa*gxa))-1;
  
  %% Efficiency Conditions
  euler   = (ct^(-sigma)) - beta*((qa*(1-delta)+ra)/qt)*(ca^(-sigma));
  qprice  = qt  - beta*gamma*(qa*(1-delta)+ra)*(1+rhoa);
  capital = Ka  - (1-delta)*Kt - ehta*(1-myu*cdft)*it;
  zcapit  = Za  - ((1+rhot)*nt-cet)/qt;
  tech    = va  - phi*vt;
  
  optc = [ euler ; qprice  ; capital ; zcapit ; tech  ];
  
  %% ----------------------%%
    xr   = [vt  Kt  ct  qt  Zt];
    REC  = [it ; Yt ; nt  ; cet]; 
    RET  = jacobian(REC,xr);
  %% ---------------------- %%
%% =========================== << Linearization >> ============================= %%   
  % Differentiation 
     xx = [ct  qt  Kt   Zt  vt   ca  qa   Ka  Za   va   ] ; % define vars
     jopt = jacobian(optc,xx);
  % Evaluate each derivative
     qt = qstar;    qa = qstar;  ct = cstar;    ca = cstar;
     Kt = kstar;    Ka = kstar;  Zt = Zstar;    Za = Zstar;
     vt = 1;        va = 1;      
   % Define Linear Coefficients  
   %  [c q K Z v]'
     coef   = double(eval(jopt));
   % In terms of % deviations from ss    
     vo = double([ cstar  qstar   kstar   Zstar   1 ]);
     TW = [ vo ;  vo  ;  vo  ;  vo ; vo  ];
   C =  [ -coef(:,1:5)   ].*TW
   B =  [ coef(:,6:10)  ].*TW
 % C[c(t)  q(t)  K(t)  Z(t)   v(t)] = B[c(t+1)  q(t+1)  k(t+1)  z(t+1)  v(t+1)]
   A = inv(C)*B   
   
  %% ---------------------------- %% 
     Recoef = double(eval(RET)); 
     ve  = double([1 kstar cstar  qstar  Zstar]);
     TQ  = [ve ; ve ; ve ;ve]
     REJ = Recoef.*TQ;
%% =============================== <<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 + e
   S1 = [ 0  0  0.1]' ;
   
   % Random diturbances
   sig = 0 ; % stop stochastic sim if sig = 0
   e = zeros(t,k)  ;
   
    
% [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')';
   
      
% [6] DRAWING FIGURES
   %% Recovering Jump variables
   %% SY = [Kt Zt vt]
   %% X  = [qt ct]
   
   XV = [X SY] ;
   qs = XV(:,1);
   cs = XV(:,2);
   Ks = XV(:,3);
   Zs = XV(:,4);
   vs = XV(:,5);
   
   XF = [vs Ks cs  qs  Zs];
   
   iYnCe = (REJ*XF')'; % i Y n ce
   XZ = [iYnCe  X];
  
  figure(1)
     subplot(3,1,1);
       plot([ XZ(:,2)  ehta.*XZ(:,1)]);
       ylabel('Y I')
       xlabel('time');
     subplot(3, 1,2);
       plot([XZ(:,3)  ehta.*XZ(:,4)]);
       ylabel('n ce')
       xlabel('time');
     subplot(3,1,3);
       plot([ cs  qs]);
       ylabel('c q')
       xlabel('time');

