% 一橋大学「金融経済論I」（２００７年度、担当：塩路)
%
% ファイル名＝islm_hyb.m
% 2007年12月1日
%
% 目的：ハイブリッド型New IS-LMモデルの分析
% モデルは教科書p.21のモデル
% 中央銀行は金利ルールをとる
% 金融緩和ショックを考える

% 背景：このプログラムのもとは加藤涼氏(IMF)のnislmq.mである。
% しかし、情報処理センターにはsymbolic math toolboxがない。
% このモデルはもともと線形なので､近似は不要である。
% よって、線形近似する部分を削除して､toolboxなしで動くようにした。
% なお、nislmq.mでは技術ショックを扱っていた。

%%  以下､加藤氏の元のコードの注釈
%%  Hybrid New ISLM model %% 
%%  By Ryo Kato, August, 2005 modified Dec. 2006

%%  This code is written by Ryo Kato (genmac@hotmail.co.jp)
%%  You can use/rewrite this code without my permission.
%%  I will not be responsible for any damage/liability 
%%  induced by running this conde.

clear all;

%% ------------------- [1] Parameter proc  ------------------------
  sigma = 1.5;   % CRRA
  alpha = 0.1;   % Slope of Phillips Curve
  q1 = 0.5;   % Taylor param 1 
  q2 = 1.5;   % Taylor param 2
  phi = 0.8;  % AR(1) in monetary policy disturbance （もとのファイルでは0.6)
  %以下がハイブリッド・パラメーター
  rho=0.1;  %加藤氏デフォルト0.91
  % Forward-looking portion param in IS (ハイブリッドでないときは１）
  beta = 0.1 ; %加藤氏デフォルト0.65
  % Forward-looking portion param in AS  (ハイブリッドでないときは１）
  
%% --------------------------[2] MODEL proc-----------------------------%%  
% このあたりは本当は不要だが､モデルの構造に関するメモとして残した。
% モデルは10変数､5本の式からなる。
% 変数は、y(t+1),pi(t+1),y(t),pi(t),v(t)・・・及び､これらの1期先の値

% Define endogenous vars  
       %('a' denotes t+1 values,'2a" denotes t+2 values.)
       %(also, '2t' is the same thing as 'a')
       %syms  y2a  pi2a  ya  pia  va  y2t pi2t yt pit vt  ; 
   % Monetary policy rule
      %ia = q1*ya+q2*pia+zetaa; 
      %ia = q1*(ya-yt)+q2*pia;  %%option speed limit policy 
  % Model equations
       %IS   = rho*y2a+(1-rho)*yt-sigma*(ia-pi2a)-ya; 
       %AS   = beta*pi2a+(1-beta)*pit+alpha*ya-pia ;
       %dum1 = ya-y2t;
       %dum2 = pia-pi2t;
       %MPs  = phi*zetat-zetaa;        
       
       %optcon  = [IS ; AS ; dum1 ; dum2; MPs];    
       
%% ------------------  [3] Linearization proc  ------------------------%%     
% 教科書p.21のモデルをそのまま入力する。
% 左辺の行列がB、右辺の行列がC
    B =  zeros(5,5) ;
    C =  zeros(5,5) ;
    B(1,1)=rho ; B(1,2)=sigma ; B(1,3)=-(1+sigma*q1) ; B(1,4)=-sigma*q2 ;
    B(2,2)=beta ; B(2,3)=alpha ; B(2,4)=-1 ; 
    B(3,3)=1 ; 
    B(4,4)=1 ;
    B(5,5)=1 ;
    C(1,3)=rho-1 ; C(1,5)=sigma*phi ;
    C(2,4)=beta-1 ;
    C(3,1)=1 ;
    C(4,2)=1 ;
    C(5,5)=phi ;
    
    % B[y(t+2)  pi(t+2)  y(t+1)  pi(t+1) zeta(t+1)] = C[y(t+1)  pi(t+1)  y(t)  pi(t) zeta(t)] 
    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 

   t = 24;			% Time span  
   % Initial Values 
   S1 = [ 0   0  -0.3 ]' ; %金融緩和ショック
   %注意：ショックはi(t+1)に与えられている。
   %i(t+1)はt期からt+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]  ; 
      X = (real(P)*SY')';
      
   % Re-definition   
     ys1 = X(:,1);  %y(t+1)の反応：こちらが見たい(ショックが起きるのはt+1期)
     pis1 = X(:,2); %pi(t+1)の反応：こちらが見たい
     ys = SY(:,1);  %y(t)の反応
     pis = SY(:,2); %pi(t)の反応
       
% [6] DRAWING FIGURES
   figure(1)
     subplot(2,1,1);
       plot( ys1 );
       ylabel('Yt')
       xlabel('time');
     subplot(2, 1,2);
       plot(pis1 );
       ylabel('dp')
       xlabel('time');

