%%  一橋大学「金融経済論I」2010年度、担当：塩路悦朗
%%  ドーンブッシュのOvershooting Model %% 
%%  ファイル名：dorn2.m
%%  作成：塩路悦朗(加藤涼氏のRBC1.mに基づく)
%%  2010年11月2日午前0時30分

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

clear all;

%% ------------------- [1] Parameter proc  ------------------------
% 変数の数を指定
njump=1 ;   %ジャンプ変数の数(為替レートのみ)
npredet=2 ; %先決変数の数(物価と貨幣)
nvar=njump+npredet ; %変数の総数

% 以下のパラメーター値を適宜設定。
a=0.01;        % 貨幣需要の利子弾力性（LMの傾きを決める)
b=0.1 ;        % 総需要の利子弾力性（ISの傾きを決める）
c=0.1 ;        % 総需要の実質為替レートに対する弾力性
d=0.01 ;       % フィリップス曲線の傾き（インフレ率がGDPギャップから受ける影響の強さ)
rho = 0.8;     % 貨幣供給のAR(1)係数

a11=d*(c+b*(1-c)/(a+b)) ;
a12=a*c*d/(a+b) ;
a21=(1-c)/(a+b) ;
a22=c/(a+b) ;

b1 = d*b/(a+b) ; 
b2 = -1/(a+b)  ; 

%% --------------------- [2] Steady State proc >> -----------------------
% 定常状態の値を計算する。(モデルが線形なので結果には無関係→省略)
  
%% --------------------------[2] MODEL proc-----------------------------%%  
   % 内生変数のリストは以下のとおり  ('a' は t+1期の変数を表す)
   % t+1期のジャンプ変数、t+1期の先決変数、t期のジャンプ変数、t期の先決変数の順
       %  ea pa ma et pt mt ; 
       %  注意:eとpは定常値からの乖離
   % 均衡条件は次のとおり
       % ea - a21*pt - (1+a22)*et -b2*mt = 0 ;
       % pa - (1-a11)*pt -a12*et -b1*mt = 0 ;
       % ma - rho*mt = 0 ;
      
%% ------------------  [3] Linearization proc  ------------------------%%     
   % 係数行列
     coef=zeros(nvar,nvar*2) ;  %係数行列は3かける6行列
     coef(1,1)=1 ; coef(1,4)=-(1+a22) ; coef(1,5)=-a21 ; coef(1,6)=-b2 ; 
     coef(2,2)=1 ; coef(2,4)=-a12 ; coef(2,5)=-(1-a11) ; coef(2,6)=-b1 ; 
     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   
     pi = SY(:,1);
     mi = SY(:,2);
     ei = X(:,1);
      
% DRAWING FIGURES
   set(0,'DefaultLineLineWidth',1.5) ; %図の線の太さを変える
   figure(1) 
       plot(0:t,mi) ;
       ylabel('変数')
       xlabel('期間');
       title('貨幣のインパルス応答関数','FontSize',14) ;
   figure(2)
     subplot(2,1,1);
       plot(0:t,pi);
       ylabel('変数')
       xlabel('期間');
       title('物価のインパルス応答関数','FontSize',14) ;
     subplot(2, 1,2);
       plot(0:t,ei);
       ylabel('変数')
       xlabel('期間');
       title('為替レートのインパルス応答関数','FontSize',14) ;