% Present-Value Test in a Markov-Switching Framework
% 
% This code implements the test implied by the present-value constraint for 
%   government debt in a Markov-switching framework. 
%   See Davig, Journal of Applied Econometrics (forthcoming) for details.
%
% Note that a minimization routine is needed (code is setup for eiter fminunc or
%   fminsearch from the Matlab optimization toolbox)
%
% Coding by :  Troy Davig
%              College of William & Mary
%              email : tadav3@wm.edu
%---------------------------------------------------------------------------- 

 
clear;
global y

% The options set parameters requited for the optimization routines used below
options = optimset('MaxFunEvals',10000000,'MaxIter',2000000000000,'TolCon',1e-12,'TolFun',1e-12);


% Read in data
%---------------------------------------------------------------------------- 
data;
y = disc_B;
[obs vars] = size(y);
%----------------------------------------------------------------------------


% Initialize
%-------------------------------------------------------------------------  

% Initial values must be specified for the optimization routine and are
% Note : p11 = exp(p) / [1 + exp(p)]
%        p11 = exp(q) / [1 + exp(q)]
%        [B*(t)|S(t)=1] = gam1 + rho B*(t-1) + e(t)
%        [B*(t)|S(t)=2] = gam2 + rho B*(t-1) + e(t)


p = 2.3;    
q = 1.8;    
gam1 = 19377;
gam2 = 73274;   
rho = 0.91;
sig = 455458311;
   
parm_guess = [p,q,gam1,gam2,rho,sig]; 
parm_guess_res = [p,q,gam2,rho,sig];    % initialization for restricted likelihood
%----------------------------------------------------------------------------


% Estimation (of unrestricted and restricted likelihoods)
%-------------------------------------------------------------------------
lik_est = ms_likfun(parm_guess,1) ;
if isnan(lik_est) == 1
    disp('Choose different initialization values...')
    break;
end
disp('maximizing (unrestricted) likelihood...')
parm_est = fminsearch('ms_likfun',parm_guess,options,1);
%parm_est = fminunc('ms_likfun',parm_guess,options,1);
lik_est = ms_likfun(parm_est,1)


lik_est_res = ms_likfun_restricted(parm_guess_res,1);
if isnan(lik_est_res) == 1
    disp('Choose different initialization values...')
    break;
end
disp('maximizing (restricted) likelihood...')
parm_est_res = fminsearch('ms_likfun_restricted',parm_guess_res,options,1);
%parm_est_res = fminunc('ms_likfun_restricted',parm_est_res,options,1);
lik_est_res = ms_likfun_restricted(parm_est_res,1) 
%-------------------------------------------------------------------------


% Compute filtered probabilities
%-------------------------------------------------------------------------
p11 = exp(parm_est(1,1)) / (1 + exp(parm_est(1,1)));
p22 = exp(parm_est(1,2)) / (1 + exp(parm_est(1,2)));
ss_prob1 = (1 - p22) / (2 - p11 - p22);
ss_prob2 = (1 - p11) / (2 - p11 - p22);

filt_probs = ms_likfun(parm_est,3);
[vals c] = size(filt_probs);
p1 = filt_probs(1:vals,1);
p2 = filt_probs(1:vals,2);
%-------------------------------------------------------------------------


% Standard Errors
%-------------------------------------------------------------------------
parn = length(parm_est);
parn_res = length(parm_est_res);
h = .00001;

for i = 1:parn;
   pert_parmest = parm_est;
   pert_parmest(1,i) = parm_est(1,i) + h;
   num_grad_i(i,:) = (ms_likfun(pert_parmest,2) - ms_likfun(parm_est,2)) / h;
end;

cov_mat = inv(num_grad_i * num_grad_i');
%-------------------------------------------------------------------------


% Output
%-------------------------------------------------------------------------
heading = ['var    ', '   estimate    ', '  ', '       se       ', '  ', 'tstat  '];
labels = ['p    : '; 'q    : '; 'gam1 : '; 'gam2 : '; 'rho  : '; 'sig  : '];
space = ['  '; '  '; '  '; '  '; '  '; '  ';];
unrestricted_estimates1 = [labels num2str(parm_est') space num2str(sqrt(diag(cov_mat))) space num2str(parm_est'./sqrt(diag(cov_mat)))];
unrestricted_estimates = [heading;unrestricted_estimates1]
p11 = exp(parm_est(1,1)) / (1 + exp(parm_est(1,1)));
p22 = exp(parm_est(1,2)) / (1 + exp(parm_est(1,2)));
trans_matrix = [p11 1-p11; p22 1-p22]

i = 1:obs;
init_yr = 1960;
t = init_yr:init_yr+obs - 1;
subplot(2,1,1)
plot(t,p2); title('Pr[S(t)=2|psi(t)]');
subplot(2,1,2)
plot(t,y); title('y(t)');
%-------------------------------------------------------------------------

chi_sq_tst = 2 * [lik_est_res - lik_est]