%% Discrete-Time Dornbusch Overshooting Model with Rational Expectations
%
%  Characteristic equation
%
% Ben J. Heijdra, Foundations of Modern Macroeconomics, Third Ed.
% Oxford University Press
%
% Chapter 5
%
% Ben J. Heijdra
% Groningen, June 9, 2015


clear
clf
clc
close all

% Structural parameters
eta         =  0.1;
delta       =  0.8;
lambda      =  0.5;
phi         =  0.2;
y_bar       =  1;
Rstar_bar   =  5;
pstar_bar   =  1;
m_bar       =  1;


% Deterministic steady-state values
y0        = y_bar
R0        = Rstar_bar
p0        = 1
m0        = m_bar;

% Euro/dollar spot rate on June 9, 2015 was 0.8869 which is 
% approximately equal to E = exp(-0.12)
e0        = -0.12; 

% Calibration parameters
alpha     =  y_bar + eta * Rstar_bar - delta *(pstar_bar + e0 - p0) ;
beta      =  m0 - p0 + lambda*Rstar_bar - y_bar ; 


% Jacobian matrix 
delta_11    = 1 + delta /(lambda+eta);
delta_12    = (1-delta) /(lambda+eta);
delta_21    = delta * lambda * phi /(lambda+eta);
delta_22    = 1 - phi * (delta*lambda + eta) /(lambda+eta);
Delta       = [delta_11 delta_12 ; delta_21 delta_22 ]
DeltaStar   = Delta - eye(1,1);


% Characteristic roots
Roots       = eig(Delta);
xi_1        = min(Roots);
xi_2        = max(Roots);

TrDelta     = delta_11 + delta_22;
DetDelta    = delta_11*delta_22 - delta_12 * delta_21 ;

s           = linspace(-1.5,3,451);
PSIs        = s.*s - s * TrDelta + DetDelta;

% Alternative phi values
% Case 2: knife-edge value for phi
phibar_L    = (lambda+eta+delta) / (delta*lambda+eta+delta);
TrDelta_L   = 2 + (delta - phibar_L * (delta*lambda + eta)) /(lambda+eta);
DetDelta_L  = (1 + delta /(lambda+eta)) * (1 - phibar_L * (delta*lambda + eta) /(lambda+eta)) ...
    - ((1-delta) /(lambda+eta)) * (delta * lambda * phibar_L /(lambda+eta)) ;
PSIs_L      = s.*s - s * TrDelta_L + DetDelta_L;

% Case 4: maximum value for phi
phibar_H    = (4*(lambda+eta)+2*delta) / (delta + 2* (delta *lambda+eta));
TrDelta_H   = 2 + (delta - phibar_H * (delta*lambda + eta)) /(lambda+eta);
DetDelta_H  = (1 + delta /(lambda+eta)) * (1 - phibar_H * (delta*lambda + eta) /(lambda+eta)) ...
    - ((1-delta) /(lambda+eta)) * (delta * lambda * phibar_H /(lambda+eta)) ;
PSIs_H      = s.*s - s * TrDelta_H + DetDelta_H;

% Case 3: high value for phi
phibar_M = phibar_H - 0.5;
TrDelta_M   = 2 + (delta - phibar_M * (delta*lambda + eta)) /(lambda+eta);
DetDelta_M  = (1 + delta /(lambda+eta)) * (1 - phibar_M * (delta*lambda + eta) /(lambda+eta)) ...
    - ((1-delta) /(lambda+eta)) * (delta * lambda * phibar_M /(lambda+eta)) ;
PSIs_M      = s.*s - s * TrDelta_M + DetDelta_M;

figure(1)
hold on;
plot(s,PSIs,'-','LineWidth',3)
box off;
hline = refline([0 0]);
plot(s,PSIs_L,'--','LineWidth',3)
plot(s,PSIs_H,':','LineWidth',3)
plot(s,PSIs_M,'-.','LineWidth',3)
plot([0,0],[-4,8],'-','LineWidth',1)
plot([1,1],[-4,8],'-','LineWidth',1)
plot([-1,-1],[-4,8],'-','LineWidth',1)
ylim([-4 7])
%legend('b=2','b=0','b=1','b=3','Location','Best')
hold off;
xlabel('s')
ylabel('\Psi (s)')
print('-deps2','-loose','F05_08')
pause(1)


% Policy functions

% Solutions for pi_p:
pi_ep_1  = (delta_11 - delta_22 - sqrt((delta_11 - delta_22)^2 + 4 * delta_12 * delta_21)) ...
    / ( 2 * delta_21)
pi_ep_2 = (delta_11 - delta_22 + sqrt((delta_11 - delta_22)^2 + 4 * delta_12 * delta_21)) ...
    / ( 2 * delta_21)
% Correct solution:
pi_ep   = min(pi_ep_1,pi_ep_2)

% Other parameters:
pi_eu   = (lambda+eta) / (lambda+eta+delta*(1-lambda*phi*pi_ep))
pi_ev   = - (delta*((1-lambda*phi*pi_ep))) / (lambda+eta+delta*(1-lambda*phi*pi_ep))
pi_e0   =  ((lambda+eta) / (delta*(1-lambda*phi*pi_ep))) ...
* ( pi_ep *(phi*(delta*lambda*pstar_bar + eta*m0) / (lambda+eta) - phi*y_bar) ...
    + Rstar_bar - (delta*pstar_bar - m0) / (lambda+eta) )

% Dynare reports the policy functions as:
%
% e_t = PI_0 + pi_p * (p_t - p_steady-state) + pi_u * u_t + pi_v * v_t 
%
% relationship with my pi_0: PI_0 = pi_0 + pi_p * P_steady-state

PI_0  =  pi_e0 + pi_ep 

pi_pp = (delta_21*pi_ep) + delta_22
pi_pu = delta_21*pi_eu
pi_pv = delta_21*(1 + pi_ev)
pi_0p = delta_21 * pi_e0 + phi*(delta*lambda*pstar_bar + eta*m0) / (lambda+eta) - phi*y_bar 
PI_0p = pi_0p + pi_pp
