%%  An New ISLM with money shocks %% 
%%  By Ryo Kato & Shinichi Nishiyama May, 2003 
%%  This code is written by Ryo Kato (kato.13@osu.edu)
%%  and Shinichi Nishiyama 
%%  You can use/rewrite this code without my permission.
%%  I will not be responsible for any damage/liability 
%%  induced by running this conde.
% Model
% IS: x(t)=x(t+1)-sigma*(i(t)-(p(t


clear all;

%% ------------------- [1] Parameter proc  ------------------------
  sigma = 1.5;   % CRRA
  kappa = 0.5;   % Cob-Dag
  beta = 0.99;   
  a = 0.5 ; 
  b = 1.4 ;   
  c = 0.3;
  rho = 0.8;


%% --------------------- [2] Steady State proc >> -----------------------
 % SS capital & ss labor
 xstar = 0;
 pstar = 0;
 istar = 0;
 zstar = 0;
 mstar = 0;
 
  
%% --------------------------[2] MODEL proc-----------------------------%%  
   % Define endogenous vars  ('a' denotes t+1 values)
       syms  xt   pt  zt  mt xa   pa  za  ma ; 
   % Eliminate LM   
       it = -(a/c)*xt+(b/c)*mt-(b/c)*pt;
   % Optimal Conditions  & state transition
       IS   = xa - sigma*it+sigma*pa - sigma*pt - xt ; 
       AS1  = beta*pa-(1+beta)*pt+zt+kappa*xt;
       AS2  = za - pt;
       Rule = rho*mt - ma;
       
       optcon  = [IS ; AS1  ; AS2 ; Rule];           
%% ------------------  [3] Linearization proc  ------------------------%%     
  % Differentiation 
     xx = [xa  pa za ma xt pt zt mt] ; % define vars
     jopt = jacobian(optcon,xx);

   % Evaluate each derivative
     xt = xstar;     xa=xstar;
     pt = pstar;     pa = pstar;
     zt = zstar;     za = zstar;
     mt = mstar;     ma = mstar;
     
   % Define Linear Coefficients  
     coef = eval(jopt);
   % In terms of % deviations from ss    
   % vo = [ xstar   pistar  vstar  ] ;
   % TW = [ vo ; vo  ; vo ] ;
 
    B =  [ -coef(:,1:4)  ] ;
    C =  [  coef(:,5:8)  ] ;
 % 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  ----------------- %%    
% [4] TIME&INITIAL VALUES 

     t = 24;			% Time span   
    
   % Initial Values 
   % state var + e
   S1 = [0  0.1]' ;

    
% [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')';
      i = -(a/c)*X(:,1)+(b/c)*SY(:,2)-(b/c)*SY(:,1);
      
% [6] DRAWING FIGURES
   figure(1)
     subplot(3,1,1);
     plot( X );
     legend('X(t)','P(t+1)',1)
       ylabel('X(t), P(t+1)')
       xlabel('time');
       subplot(3, 1, 2);
       plot(i);
       ylabel('i');
       xlabel('time');
     subplot(3, 1, 3);
     plot( SY );
     legend('P(t)','m(t)',1)
       ylabel('State')
       xlabel('time');
       
 



