% 一橋大学「金融経済論I」（２００７年度、担当：塩路)
%
% ファイル名＝islm_rho.m
% 2007年12月1日
%
% 目的：New IS-LMモデルの分析
% islm_int.mの拡張
% Calvo parameterであるρが変わったら、金融政策に対する反応は？
% 中央銀行は金利ルール(テイラールール)をとる

% 背景：このプログラムのもとは加藤涼氏(IMF)のnislm.mとnislmq.mである。

clear all;

%% ------------------- [0] 準備  ------------------------
%rhovec=0.1:0.1:0.9 ; %rhoの値を色々変えてみる
rhovec=[0.1 0.5 0.9] ;
t = 24;			     % 何期先の反応まで見たいか(ショックのあった期を0期とする)

xsize=size(rhovec) ;
xnum=xsize(2) ;
yvec=zeros(t+1,xnum) ; %結果を保存しておく入れ物を作っておく
pivec=yvec ;
ivec=yvec ;
mvec=yvec ;

%% ------------------- [1] Parameter proc  ------------------------
for ix=1:xnum ; %ここからループが始まる

  sigma = 1.5;   % CRRA、教科書では1/θ
  theta=1/sigma ;
  lambda=1 ;
  myu = 2;       % preference
  phi = 0.8;     % AR(1) in Taylor rule error term
  q1 = 0.5 ;     % Taylor param 1 
  q2 = 1.5 ;     % Taylor param 2
 
  rho=rhovec(ix) ;
  alpha = rho^2*(lambda+theta)/(1-rho);   % Slope of Phillips Curve、75ページ

%% ------------------  [3] Linearization proc  ------------------------%%     
% 教科書p.75-76の記述からモデルを求めそのまま入力。
% 左辺の行列がB、右辺の行列がC
    B =  zeros(3,3) ;
    C =  zeros(3,3) ;
    B(1,1)=1 ; B(1,2)=sigma ; 
    B(2,2)=1 ; 
    B(3,3)=1 ;
    C(1,1)=1+sigma*q1 ; C(1,2)=sigma*q2 ; C(1,3)=sigma*phi ;
    C(2,1)=-alpha ; C(2,2)=1 ;
    C(3,3)=phi ;
    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 
     
  %  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  ----------------- %%    
% TIME&INITIAL VALUES 

   % Initial Values 
   S1 = [ -0.3 ]' ; %先決変数はzetaだけ、教科書にならって利下げショックを考える

% 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   
     ys = X(:,1);     %GDPギャップの反応
     pis = X(:,2);    %インフレ率の反応
     zetas = SY(:,1); %金融政策を表す項ζ
     
     is = q1*ys+q2*pis+zetas ; %金利の反応をリカバー
     ms = -(1/myu)*is+(1/(myu*sigma))*ys ; %実質貨幣の反応をリカバー

   % 結果の保存
     yvec(:,ix)=ys ; pivec(:,ix)=pis ; ivec(:,ix)=is ; mvec(:,ix)=ms ;
     
end ; %ここでループが終わる     
     
     
% [6] DRAWING FIGURES
  set(0,'DefaultLineLineWidth',1.5) ; %図の線の太さを変える
  figure(1)
     subplot(2,2,1);
       plot(yvec); title('GDPギャップの反応','FontSize',14) ;  
     subplot(2,2,2);
       plot(pivec); title('インフレ率の反応','FontSize',14) ;  
     subplot(2,2,3);
       plot(ivec); title('名目利子率の反応','FontSize',14) ;  
     subplot(2,2,4);
       plot(mvec); title('実質貨幣残高の反応','FontSize',14) ;  
 