% 一橋大学「金融経済論I」（2012年度、担当：塩路)
%
% ファイル名＝ps3_2012.m
% 2012年11月28日
%
% 目的：課題3のこたえ

% このファイルを動かすためには、同じディレクトリ内に次が必要である。
%   fdjac_japan.m: 数値微分を取るためのプログラム(Fackler and Mirandaにもとづく）
%   ps3_2012_model.m: モデルは、この中で定式化されている
%   ps3_2012_others.m: モデルの内生変数からGDPを計算するプログラム

clear all;

% 変数の数を指定
njump=2 ;   %ジャンプ変数の数
npredet=2 ; %先決変数の数
nvar=njump+npredet ; %変数の総数


%% ------------------- [1] Parameter proc  ------------------------
  global sigma alpha myu beta lam psi phi ;
  sigma = 1.5;   % CRRA
  alpha = 0.3;   % Cobb-Doug
  beta = 0.99;   % discount factor
  lam=3 ; %elasticity of depreciation to capacity utilization
  deltastar=0.025 ; %steady state depreciation
  psi=(1/beta-1+deltastar)/alpha ; %importance of endogenous capacity utilization
  myu = psi/(1+lam)-deltastar; % constant term in depreciation
  phi = 0.95;     % AR(1) in tech

%% --------------------- [2] Steady State proc >> -----------------------
zstar=1 ;
vstar=1;
kstar=(alpha/(1/beta-1+deltastar))^(1/(1-alpha)) ;
Ystar=kstar^alpha ;
cstar=Ystar-deltastar*kstar ;
  
%% --------------------------[2] MODEL proc-----------------------------%%  
   % 内生変数のリスト  ('a' denotes t+1 values)
   % t+1期のジャンプ変数、t+1期の先決変数、t期のジャンプ変数、t期のジャンプ変数の順
   % このモデルではジャンプ変数はzとcで先決変数はkとv
       %  za ca   ka va  zt ct   kt vt   
    
%% ------------------  [3] Linearization proc  ------------------------%%     
% DIFFERENTIATION OF THE MODEL
% Specify the point at which the derivatives are taken
    xss=[zstar cstar kstar vstar zstar cstar kstar vstar] ;
% Differentiation 
    coef1 = fdjac_japan('ps3_2012_model',xss);
    
% DIFFERENTIATION OF THE GDP EQUATION
% Specify the point at which the derivatives are taken
    yss=[zstar  kstar  vstar] ;
% Differentiation 
    coy1 = fdjac_japan('ps3_2012_others',yss);

% TRANSFORMATION INTO A LOG LINEARIZED MODEL
% In terms of % deviations from ss    
    vo = [ zstar cstar   kstar   vstar  ] ;
    TW = [ vo ; vo  ; vo ; vo ] ;
 
    B =  [ -coef1(:,1:4)  ].*TW ;
    C =  [  coef1(:,5:8)  ].*TW ;
        % B[z(t+1)  c(t+1)  k(t+1)  v(t+1)] = C[z(t)  c(t)  k(t)  v(t)]
    A = inv(C)*B  %(Linearized reduced form )
   
   % For GDP( optional) 
      ve  = [zstar  kstar  vstar  ];
      NOM = [Ystar  Ystar  Ystar  ];
      PPX = coy1.*ve./NOM;

%% =========== [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  ----------------- %%    
% モデルが一義的な解を持つかどうかをチェック(Blanchard-Khan条件)
 if njump~=n ;
     disp('ジャンプ変数の数と発散的な固有値の数が異なります。') ;
     disp('計算を停止します') ; return ;
 end ;

% [4] TIME&INITIAL VALUES 

     t = 24;			% Time span   
    
   % Initial Values 
   % state var + e
   S1 = [ 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')';
      
   % Re-definition   
     zi = X(:,1);
     ci = X(:,2);
     ki = SY(:,1);
     vi = SY(:,2);
      
     XI = [zi ki vi];
     Yi = (PPX*XI')';      
% [6] DRAWING FIGURES
   set(0,'DefaultLineLineWidth',1.5) ; %図の線の太さを変える
   figure(1)
     subplot(2,1,1);
       plot( Yi );
       ylabel('Output')
       xlabel('time');
     subplot(2, 1,2);
       plot([ci zi] );
       ylabel('Consumption(bl)  Utilization(grn)')
       xlabel('time');