%% Kiyotaki-Moore Model (Credit Cycles)
%%  JPE (1997)
%%  Ryo Kato September 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.


clear all;
close all;

%% -----------------------------------------------------------
% [1]  Parameters
  a   = 1;
  phi = 20;
  R   = 1.01;
  lamda = 0.975;
  pai = 0.1 ;  
  v  = 2;
  
%% ============================= << SS Values >> ================================ %%  

 DKstar = (a + lamda*phi - phi)/(R-1);
 qstar  = (R/(R-1))*(pai*a - (1-lamda)*(1- R + pai*R)*phi )/(lamda*pai + (1-lamda)*(1- R + pai*R));
 xstar  = qstar * (R-1)/R;
 Kstar  = xstar + v ;
 Dstar  = DKstar * Kstar;
  
 qKD = [qstar  Kstar  Dstar]'
                     
%% ============================= << MODEL >> ================================ %%  
 
 syms  Dt Kt qt Da  Ka qa  ;  
                                              
  capital =  (pai/(qt + phi -qa/R))*( (a+lamda*phi+qt)*Kt - R*Dt) + (1-pai)*lamda*Kt - Ka ;
  debt    =  R*Dt + qt*(Ka - Kt) + phi*(Ka-lamda*Kt) - a*Kt -Da ;
  euler   =  qt - qa/R - (Kt - v) ;
  
  optc = [ euler ;  debt ; capital  ];
  
%% =========================== << Linearization >> ============================= %%   
  % Differentiation 
     xx = [qt  Dt  Kt   qa  Da   Ka   ] ; % define vars
     jopt = jacobian(optc,xx);
  % Evaluate each derivative
     qt = qstar;    qa = qstar;  
     Dt = Dstar;    Da = Dstar;
     Kt = Kstar;    Ka = Kstar;  
   % Define Linear Coefficients  
     coef   = double(eval(jopt));
   % In terms of % deviations from ss    
     vo = double([ qstar  Dstar   Kstar  ]);
     TW = [ vo ;  vo  ;  vo   ];
   C =  [ -coef(:,1:3)   ].*TW;
   B =  [ coef(:,4:6)  ].*TW;
 % C[q(t)  D(t)  K(t)] = B[q(t+1)  D(t+1)  K(t+1)]
   A = inv(C)*B   
   
  
%% =============================== <<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 = 80;			% Time span you want to simulate   
    
   % Initial Values 
   % state var + e
   S1 = [ 0  0.01]' ;
   
   % 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 = [Dt Kt]
   %% X  = [qt]
   
   qs = X(:,1);
   Ks = SY(:,2);
   Ds = SY(:,1);
   
 
  figure(1)
     subplot(2,1,1);
       plot( qs  );
       ylabel('qt')
       xlabel('time');
     subplot(2, 1,2);
       plot([Ks  Ds]);
       ylabel('Kt Dt')
       xlabel('time');



