%%  一橋大学「金融経済論I」2012年度、担当：塩路悦朗
%%  トービンのqモデル %% 
%%  ファイル名：tobinq2.m
%%  作成：塩路悦朗(加藤涼氏のRBC1.mに基づく)
%%  2012年11月6日午後10時50分

%%  目的：講義資料にあるトービンのqモデルに、Bの項がAR(1)過程に従うという
%%  仮定を追加した。
%%  このコードは、第0期にBに1単位のショックが加わったとして、
%%  その後内生変数（資本ストックとq)がどのように変化するかを分析する。
%%  （内生変数のショックに対する反応を通時的に追ったものをインパルス応答関数という)

clear all;

%% ------------------- [1] Parameter proc  ------------------------
% 変数の数を指定
njump=1 ;   %ジャンプ変数の数(qのみ)
npredet=2 ; %先決変数の数(資本ストックとB)
nvar=njump+npredet ; %変数の総数

% 以下のパラメーター値を適宜設定。
% 以下のパラメーター値を適宜設定。
A=1;           
r=0.05 ;
psi=10 ;
rho = 0.95;     % BのAR(1)係数

%% --------------------- [2] Steady State proc >> -----------------------
% 定常状態の値を計算する。(モデルが線形なので結果には無関係→省略)
  
%% --------------------------[2] MODEL proc-----------------------------%%  
   % 内生変数のリストは以下のとおり  ('a' は t+1期の変数を表す)
   % t+1期のジャンプ変数、t+1期の先決変数、t期のジャンプ変数、t期の先決変数の順
       %  qa ka ba qt kt bt ; 
       %  注意:qとkは定常値からの乖離
   % 均衡条件は次のとおり
       % -qa+A*ka=-(1+r)*qt+bt ;
       % ka=(1/psi)*qt+kt ;
       % ba - rho*bt = 0 ;
      
%% ------------------  [3] Linearization proc  ------------------------%%     
   % 係数行列
     coef=zeros(nvar,nvar*2) ;  %係数行列は3かける6行列
     coef(1,1)=-1 ; coef(1,2)=A ; coef(1,4)=1+r ; coef(1,6)=-1 ;  
     coef(2,2)=1 ; coef(2,4)=-1/psi ; coef(2,5)=-1 ;  
     coef(3,3)=1 ; coef(3,6)=-rho ; 
     
    B =  [  coef(:,1:nvar)  ] ; % t+1期の変数の係数
    C =  [ -coef(:,nvar+1:2*nvar)  ] ; % t期の変数の係数
    % この時点で、モデルはB*x(t+1)-C*x(t)=0という形をしている。
    A = inv(C)*B  % ここでA*x(t+1)=x(t)という形のモデルに直した。
   
%% =========== [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 
     
  % 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 ;
     
% TIME&INITIAL VALUES 

     t = 200;			% 何期後の反応まで求めたいか
    
   % 先決変数の初期値を与える＝モデルにショックを加える 
       S1 = [ 0 1 ]' ;
   % 第2の先決変数であるmに1を加える＝1単位の正の貨幣供給ショック

    
% SIMULATION 
      Ss = S1;
      S = zeros(t,k) ;  %以下で先決変数の反応を求める
        for i = 1:t
          q = AA*Ss  ;
   	      S(i,:) = q';
          Ss = S(i,:)';
     	end
      SY = [S1' ;S]  ; %SYが先決変数の反応
      X = (real(P)*SY')'; %Xがジャンプ変数の反応
      
   % Re-definition   
     ki = SY(:,1);
     bi = SY(:,2);
     qi = X(:,1);
      
% DRAWING FIGURES
   set(0,'DefaultLineLineWidth',1.5) ; %図の線の太さを変える
   figure(1) 
       plot(0:t,bi) ;
       ylabel('変数')
       xlabel('期間');
       title('Bのインパルス応答関数','FontSize',14) ;
   figure(2)
     subplot(2,1,1);
       plot(0:t,ki);
       ylabel('変数')
       xlabel('期間');
       title('資本ストックのインパルス応答関数','FontSize',14) ;
     subplot(2, 1,2);
       plot(0:t,qi);
       ylabel('変数')
       xlabel('期間');
       title('qのインパルス応答関数','FontSize',14) ;