% 一橋大学「金融経済論I」（２０１０年度、担当：塩路)
%
% ファイル名＝rbc1.m
% 2010年11月9日
% 授業後訂正あり
%
% 目的：RBCモデルの分析を、手で対数線形近似をした上で、行う

% 背景：このプログラムの基本は加藤涼氏(IMF)のrbc1.mである。
% しかし、加藤氏のプログラムは、Matlabの別売りソフトである、
% symbolic math toolbox(数万円)なしではモデルの対数線形近似を取ってくれない。
% このソフトは我らが情報基盤センターには、ない。
% そこで、パソコンで微分を取ってもらう代わりに地道に自分の手で微分をとるようにすることで、
% このtoolboxなしで分析できるようにすることを目指すことにした。


%%  加藤氏の元のコードの注釈
%%  Basic RBC model %% 
%%  By Ryo Kato, modified January, 2004
%%               modified June 2005
%%  variable elastic labor supply
%%  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.


%%  ATTENTION!!
%%  Labor disutility ==>   -myu*(lt^(1+lamda))/(1+lamda)
%%  教科書第1章とやや定式化が異なっているので注意(第7章と同じ)

clear all;       %これまでの計算結果は忘れなさい

%% ------------------- [1] パラメーターの設定  ------------------------
   % 変数の数を指定
   njump=2 ;   %ジャンプ変数の数
   npredet=2 ; %先決変数の数
   nvar=njump+npredet ; %変数の総数

  %パラメーター値は教科書P.38に従う
  sigma = 1.5 ;  % 効用関数のパラメーター、CRRA
  alpha = 0.3;   % 資本分配率
  myu = 1        % 労働の負効用の係数
  beta = 0.99;   % 主観的割引因子(ファクター)
  delta = 0.025; % 資本減耗率
  lamda = 1 ;    % 労働供給の弾力性 
  phi = 0.8;     % 生産性のAR(1)係数、教科書には0.9と書いてあるが、どうも0.8を使ったのではないか？？

%% --------------------- [2] 非確率的定常状態の導出 >> -----------------------
 % SS capital & ss labor
 % (1) real rate  (By SS euler)
        kls = (((1/beta)+delta-1)/alpha)^(1/(alpha-1)); %資本・労働比率
 % (2) wage
        wstar = (1-alpha)*(kls)^alpha;                  %実質賃金
 % (3) Labor and goods market clear
        clstar = kls^alpha - delta*kls;                 %消費・労働比率
        lstar = ((wstar/myu)*(clstar^(-sigma)))^(1/(lamda+sigma)) ; %労働
        kstar = kls*lstar ;                             %資本ストック
        cstar = clstar*lstar ;                          %消費
        ckstar= cstar/kstar ;                           %消費・資本比率
        vstar = 1;                                      %生産性の定常値
        Ystar = (kstar^alpha)*(lstar^(1-alpha));        %生産
     
    ssCKoLY = [ cstar  kstar ; lstar  Ystar ]           %定常値の表示
   
  
%% --------------------------[3] モデルの設定-----------------------------%%  
   % 内生変数のリスト  ('a' denotes t+1 values)
   % t+1期のジャンプ変数、t+1期の先決変数、t期のジャンプ変数、t期のジャンプ変数の順
   % このモデルではジャンプ変数はlとcで先決変数はkとv

   % ここでは係数行列を手で求めて入力する。詳細は講義ノートに。
   % 講義ノートが間違っていたので、一部に訂正があります。
   % このコードそのものにも1か所誤りがありました。
   
   a=beta*(1/beta-1+delta)/sigma ; %一部訂正（講義ノートの誤り）
   b=kls^(-(1-alpha)) ;
   c=ckstar ;
   
   coef1left=[-(1-alpha)*a 1 (1-alpha)*a -a ; 0 0 0 0 ; 0 0 1 0 ; 0 0 0 1] ; %一部訂正（講義ノートの誤り）
   coef1right=[0 -1 0 0 ; -(alpha+lamda) -sigma alpha 1 ; -(1-alpha)*b c -(1-delta+alpha*b) -b ; 0 0 0 -phi] ; %一部訂正（講義ノートの誤り＋コード自体の誤り）
   
   coef1=[coef1left coef1right] ;
   
   B =  [ -coef1(:,1:nvar)  ] ;
   C =  [  coef1(:,nvar+1:2*nvar)  ] ;
        % B[l(t+1)  c(t+1)  k(t+1)  v(t+1)] = C[l(t)  c(t)  k(t)  v(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 
     
  % [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]  シミュレーションを行う  ----------------- %%    
% [4] シミュレーション期間を設定し初期値を設定する 

    % モデルが一義的な解を持つかどうかをチェック(Blanchard-Khan条件)
     if njump~=n ;
         disp('ジャンプ変数の数と発散的な固有値の数が異なります。') ;
         disp('計算を停止します') ; return ;
     end ;

    t = 24 ;			% シミュレーション期間   
    
     S1 = [ 0   1 ]';   %初期値：kは定常状態、生産性に1単位のショックが加わったとする。

    
% [5] シミュレーションを行う
      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')'; %ジャンプ変数のシミュレーション結果
      
   % 結果を取り出す
     li = X(:,1);
     ci = X(:,2);
     ki = SY(:,1);
     vi = SY(:,2);
      
% [6] 結果をグラフ化：教科書P.39に合わせる
   set(0,'DefaultLineLineWidth',1.5) ; %図の線の太さを変える
   figure(1)
     subplot(2,2,1);
       plot( [0 ; li] ); xlim([0 24]) ;
       ylabel('労働')
       xlabel('期間');
     subplot(2,2,2);
       plot( [0 ; ci] ); xlim([0 24]) ;
       ylabel('消費')
       xlabel('期間');
     subplot(2,2,3);
       plot( [0 ; vi] ); xlim([0 24]) ;
       ylabel('生産性')
       xlabel('期間');
     subplot(2,2,4);
       plot( [ 0 ; ki] ); xlim([0 24]) ;
       ylabel('資本ストック')
       xlabel('期間');
