% 一橋大学「金融経済論I」（2012年度、担当：塩路)
%
% ファイル名＝kiyotakimoore1.m
% 2013年1月22日
%
% 目的：Kiyotaki-Mooreモデルの分析
%
% 経済が生産性ショックに見舞われたとする
% ここで生産性ショックとは、パラメーターaに対するショック

%% ------------------- [1] パラメーターの設定  ------------------------
clear all;       %これまでの計算結果は忘れなさい

% 変数の数を指定
   njump=1 ;   %ジャンプ変数の数:qだけ
   npredet=3 ; %先決変数の数:K、Bと生産性a
   nvar=njump+npredet ; %変数の総数
% このファイルを動かすためには、同じディレクトリ内に次が必要である。
%   fdjac_japan.m: 数値微分を取るためのプログラム(Fackler and Mirandaにもとづく）
%   kiyotakimoore1_model.m: モデルは、この中で定式化されている

  global R v pai lamda phi astar psi  ;
 % パラメーター値の設定
  astar   = 1; %生産性パラメーターaの定常値
  phi = 20;
  R   = 1.01;
  lamda = 0.975;
  pai = 0.1 ;  
  v  = 2;
  psi=0.9 ; %ショックのAR1パラメーター

%% --------------------- [2] 定常状態の計算 -----------------------
%生産性aの値がastarに等しいときの経済の定常状態を求める。 
 BKstar = (astar + lamda*phi - phi)/(R-1);
 qstar=(pai*astar-(1-lamda)*phi-pai*R*BKstar)/((1-(1-pai)*lamda)*(1-1/R)-pai) ;
 xstar  = qstar * (R-1)/R;
 Kstar  = xstar + v ;
 Bstar  = BKstar * Kstar;
  
%% -----------------------[3] モデルの線形近似を取る -----------------------------%%  
   % 内生変数のリスト  ('a' denotes t+1 values)
   % t+1期のジャンプ変数、t+1期の先決変数、t期のジャンプ変数、t期のジャンプ変数の順
   % このモデルではジャンプ変数はq(t)で先決変数はk(t-1)とb(t-1)とa(t)
   % バージョン1ではこの順序が間違っていた。
       %  qa ka ba aa   qt kt bt at   
    
% モデルの各式を各内生変数について微分する
% どのポイントで微分するか
    xss=[qstar Kstar Bstar astar qstar Kstar Bstar astar] ;
% どの関数を微分するか
    coef1 = fdjac_japan('kiyotakimoore1_model',xss);
    
% 線形近似から対数線形近似へ
% こうすることで、式は定常状態からの％乖離の形で表現できる
    vo = [ qstar Kstar   Bstar   astar  ] ;
    TW = [ vo ; vo  ; vo ; vo ] ;
 
    B =  [ -coef1(:,1:4)  ].*TW ;
    C =  [  coef1(:,5:8)  ].*TW ;
        % B[q(t+1)  K(t)  B(t)  a(t+1)] = C[q(t)  K(t-1)  B(t-1)  a(t)]
    A = inv(C)*B  %(Linearized reduced form )
   
%% =========== [4] モデルを解く（ここは加藤涼氏のコードのまま、手をつけない) ============== %%
%  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]  モデルの解をもとにシミュレーションを行う  ----------------- %%    
% モデルが一義的な解を持つかどうかをチェック(Blanchard-Khan条件)
 if njump~=n ;
     disp('ジャンプ変数の数と発散的な固有値の数が異なります。') ;
     disp('計算を停止します') ; return ;
 end ;

 t = 99;			% 99期先までのインパルス応答関数を求める (ショックのあった期を入れると100期）  
    
% 初期値：KとBはもともと定常値にあった。急にaが定常値から1%増加した。 
   S1 = [ 0 0 1 ]' ; %ジャンプしない変数の初期値（定常値からの乖離）
    
% シミュレーションを走らせる 
      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')'; %ジャンプする変数に関するシミュレーション結果
      
   % 結果を取り出して名前を付ける
     qi = X(:,1);
     ki = SY(:,1);
     bi = SY(:,2);
     ai = SY(:,3);
      
% インパルス応答関数の図を描く
   set(0,'DefaultLineLineWidth',1.5) ; %図の線の太さを変える
   figure(1)
     subplot(2,2,1);
       plot( ai );
       title('生産性')
     subplot(2,2,2);
       plot( qi );
       title('地価')
     subplot(2,2,3);
       plot( ki );
       title('K（耕地）')
     subplot(2,2,4);
       plot( bi );
       title('負債')

