% BASIC Monetary Policy Model (Standard MODEL for OPEN Economy)
% Dec 1 2001

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

%%  This code is a beta version.    

% Persistent inflation with persistent policy, government ex and 
% Exchange rate rule  DS(t+1) = (a+1)*Pi[t+1] + b*(y[t+1]-y*[t+1])


%% ------------------- << Stochastic Similation >> --------------------- %%

% # of simulations = # of PLOTS 
lmax = 300 ;
% Bound on policy parametar 
low = 0.1 ;
upper = 2 ;

% Stochastic simulations
  arange = linspace(low,upper,lmax);
  brange = linspace(low,upper,lmax);
  RESULT = [] ; % Storing results
  stdv =[];     % Storing standard deviations of [gap pai call] 

for l = 1 : lmax   

%% ------------------- <<MODEL>>--------------------------------------%%
% [1] DEFINITION of the MODEL
% C*[y(t+1) y*(t+1) pi(t+1) R(t) s(t+1) e(t+1) v(t+1) ]'
%  =B[y(t) y*(t) pi(t) R(t-1) s(t)  e(t) v(t)]'  

% s(t) = pi(t-1)  v= gov exp  e = supply shock
lamda = 2; % Exchange rate and POTGDP
w1 = 0.6 ;   % weight on c 
w2 = 0.1 ;   % weight on export
w3 = 0.15 ;   % weight on gov exp
sigma = 0.1; % IES
beta = 0.5 ; % pi lead
a = arange(l) ;    
b = brange(lmax+1-l) ;
k = 0.02 ;    % Phil slope
phi = 0.6 ;  % supply shock persistence
rho = 0.8 ;  % policy persistence
chi = 0.6 ;  % gov exp shock persistence

% Matirx
c11	=	1;	c12	=	lamda*w2	;	c13	=	sigma*w1	;	c14	=	-sigma*w1	; c15=0 ; c16=0 ; c17=-w3 ;
c21	=	0;	c22	=	lamda	;	c23	=	-1	;	c24	= 1	;	c25	=	0	;	c26	=	0	;	c27	=	0;
c31	=	0;	c32	=	0	;	c33	=	beta	;	c34	= 0	;	c35	=	0	;	c36	=	0	;	c37	=	0;
c41	= -b;	c42	=	b  ;	c43	=	-1-a	;	c44	=	1	;	c45	=	0	;	c46	=	0	;	c47	=	0;
c51	=	0;	c52	=	0	;	c53	=	0	;	c54	=	0	;	c55	=	1	;	c56	=	0	;	c57	=	0;
c61	=	0;	c62	=	0	;	c63	=	0	;	c64	=	0	;	c65	=	0	;	c66	=	1	;	c67	=	0;
c71	=	0;	c72	=	0	;	c73	=	0	;	c74	=	0	;	c75	=	0	;	c76	=	0	;	c77	=	1;
  
  C = [c11 c12 c13 c14 c15 c16 c17 ; c21 c22 c23 c24 c25 c26  c27; c31 c32 c33 c34 c35 c36 c37 ; 
       c41 c42 c43 c44 c45 c46  c47; c51 c52 c53 c54 c55 c56 c57 ; c61 c62 c63 c64 c65 c66  c67; 
       c71 c72 c73 c74 c75 c76 c77;  ] ;
       
b11	=	1;	b12	=	lamda*w2 ;	b13	=	0;	b14	=	0	;	b15	=	0	;	b16	=	0	;	b17	=	-w3;
b21	=	0;	b22	=	lamda	 ;	b23	=	0;	b24	=	0	;	b25	=	0	;	b26	=	0	;	b27	=	0;
b31	= -k ;	b32	=	k	;	b33	=	1	;	b34	=	0	;	b35	= beta-1;   b36	= -1    ;	b37	=	0;
b41	= -b;	b42	=	b	;	b43	= -1-a	;	b44	=	1	;	b45	=	0	;	b46	=	0	;	b47	=	0;
b51	=	0;	b52	=	0	;	b53	=	1	;	b54	=	0	;	b55	=	0	;	b56	=	0	;	b57	=	0;
b61	=	0;	b62	=	0	;	b63	=	0	;	b64	=	0	;	b65	=	0	;	b66	=	phi	;	b67	=	0;
b71	=	0;	b72	=	0	;	b73	=	0	;	b74	=	0	;	b75	=	0	;	b76	=	0	;	b77	=	chi;
  
  B = [b11 b12 b13 b14 b15 b16 b17 ; b21 b22 b23 b24 b25 b26 b27 ; b31 b32 b33 b34 b35 b36 b37 ; 
       b41 b42 b43 b44 b45 b46 b47 ; b51 b52 b53 b54 b55 b56 b57 ; b61 b62 b63 b64 b65 b66 b67; 
       b71 b72 b73 b74 b75 b76 b77;  ] ;

% Reduced Form X(t+1) = A*X(t)
A = inv(C)*B;

% Reduced Form X(t+1) = A*X(t)
A = inv(C)*B;


% --------------------------------<<Solution proc>> -------------------- %
% [2] EIGEN DECOMPOSITION
   [W V] = eig(A);
   Q = inv(W);
   W*V*Q; 
   theta = diag(V);
  % Extract unstable vectors
      UQ = [];
      jw = 1;
      for j = 1:length(theta)
      	if abs(theta(j)) > 0.999
         	UQ(jw,:) = Q(j,:);
         	jw = jw+1;
      	end
      end
   % Extract stable vectors
     SQ = [];
     jjw = 1;
      for jj = 1:length(theta)
      	if abs(theta(jj)) <0.999
         	SQ(jjw,:) = Q(jj,:);
         	jjw = jjw+1;
      	end
      end
    % Extract stable roots
      VLL = [];
      jjjw = 1;
      for jjj = 1:length(theta)
      	if abs(theta(jjj)) <0.999
         	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)
      VLLL = [];
   hh = 1;
   for h = 1:k
      VLLL(:,hh) = VLL(h)*ones(k,1);
      hh = hh +1 ;
   end
   VLLL;   
   VL = eye(k).*VLLL  ;
   
   % 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)  ; 
    
    
%% --------------------------<< SIMULATION PROGRAM  >>----------------------------- %%    
% [4] TIME&INITIAL VALUES 

     t = 200;			% Time span you want to simulate   
    
   % Initial Values 
   % state var + e
   S1 = [zeros(k-2,1) ; 0 ; 0 ] ;
   
   % Random diturbances
   sig = 0.5 ; % stop stochastic sim if sig = 0
   gsig = 0.6 ;
   RNG = randn(t,1,lmax);
   RNSP = gsig*randn(t,1,lmax);
   e = sig*[zeros(t,k-2,lmax)  RNSP  RNG ] ;
   
        
% [5] SIMULATION 
      Ss = S1;
      S = zeros(t,k) ;  
        for i = 1:t
          q = AA*Ss + e(i,:,l)' ;
   	      S(i,:) = q';
          Ss = S(i,:)';
     	end
      SY = [S1' ;S]   ;
      X = (real(P)*SY')';
  % GDP gap    
      gap = X(1:t-1,1)-X(1:t-1,2);  % 
      pai = X(1:t-1,3);               %  
      call =SY(2:t,1);          % R(-1) is state
      
      RESULT(:,:,l)=[gap pai call];
      stdv(l,1:3) = [std(RESULT(:,:,l))];
      
   end

   stdv;
   avarY = stdv(:,1).^2;
   avarPI = stdv(:,2).^2;
   avarCR = stdv(:,3).^2;
   
   ALT = [avarY avarPI  avarCR]
   
 % [6] Drawing Pictures  
    figure(1)
   plot(avarY, avarPI,'r*');
 %  plot(avarY, avarPI,'r*',varY, varPI,'b+');
    xlabel('Output gap Variance');
    ylabel('Inlation Variance');
    
     
   
   
   
   
      
