/****************************************************************************/;
/****************************************************************************/;
/*****   Solve-the-Equation Plug-In (SP) Bandwidth Choice Rule for      *****/;
/*****   HAC Estimation with the Bartlett Kernel                        *****/;
/****************************************************************************/;
/*****   Source: Hirukawa, M. (2006): "A Two-Stage Plug-In Bandwidth    *****/;
/*****   Selection and Its Implementation for Covariance Estimation."   *****/;
/****************************************************************************/;
/*****   Last Modified: May 11, 2006                                    *****/;
/****************************************************************************/;
/****************************************************************************/;



new;
format 6,4;

/****************************************************************************/;
/*****   1. Defining observation matrix and weighting vector            *****/;                  
/****************************************************************************/;

@-----------------------------------------------------------------------------@
@   DEFINE "G", THE MATRIX OF ORTHGONALITY CONDITIONS.  THE PROCESS {g(t)}    @
@   IS ASSUMED TO BE A ZERO-MEAN PROCESS.                                     @
@-----------------------------------------------------------------------------@

t = 200;    @ number of observations: INPUT HERE. @
d = 2;      @ dimension of vector process {g(t)}: INPUT HERE. @
load g[t,d] = c:\iphac\gauss\test.dat;  @ t x d matrix of orthogonality conditions: INPUT HERE. @ 

@-----------------------------------------------------------------------------@
@   DEFINE "w", THE WEIGHTING VECTOR.                                         @
@-----------------------------------------------------------------------------@

w = 0|ones(d-1,1);  @ w = (0,1,...,1)' is assumed, where the first row corrensponds to the intercept: CHANGE WHENEVER NECESSARY. @ 

@-----------------------------------------------------------------------------@
@   DO NOT CHANGE THE VALUES BELOW.                                           @
@-----------------------------------------------------------------------------@
				
grid = .25;			@ step size for grid search @
lb = 0;             @ minimum possible bandwidth @ 
ub = t;			    @ maximum possible bandwidth: setting equal to the sample size usually works. @ 
maxtol = 10^(-6);   @ maximum tolerance for fixed-point iteration @

q = 1;              @ characteristic exponent of Bartlett kernel @	
kq = 1;             @ generalized derivative of Bartlett kernel at origin @
intk2 = 2/3;        @ squared integral of Bartlett kernel @
intx2k2 = 1/15;     @ "2nd-order moment" of Bartlett kernel @



/****************************************************************************/;
/*****   2. Implementing Bartlett-IP bandwidth choice rule              *****/;                  
/****************************************************************************/;

@-----------------------------------------------------------------------------@
@   OBTAIN THE PROXY OF ALPHA(Q), THE UNKNOWN QUANTITY IN IP ALGORITHM.       @
@-----------------------------------------------------------------------------@

h = g*w;        
h0 = h[2:t];
h1 = h[1:t-1];
phi = invpd(h1'h1)*(h1'h0);     @ least squares estimate of AR(1) coefficient @
alpha = (phi^2+1)/(phi^2-1);	@ proxy of unknown quantity in IP algorithm @

@-----------------------------------------------------------------------------@
@   FIXED POINT ITERATION USED HERE IS A COMBINATION OF GRID SEARCH &         @
@   BISECTION SEARCH.  FIRST, GRID SEARCH ROUTINE IS GIVEN BELOW.             @
@-----------------------------------------------------------------------------@

sthat = ub; @ estimate of optimal bandwidth @

print "Grid search starts here...";
print "";
print "Lower bound = " lb;;
print "Upper bound = " ub;
print "";

bt = ((alpha^2)*intk2*ub^(2*q+1)/((2*q+1)*intx2k2))^(1/(4*q+1));
{rn} = bartlett(h,bt,q);    @ estimate of q'th generalized derivative of spectral density at origin @
{rd} = bartlett(h,bt,0);    @ estimate of spectral density at origin @
r2 = (rn/rd)^2;   @ squared normalized curvature @
fub = ub-(q*(kq^2)*r2*t/intk2)^(1/(2*q+1));

j = 1; do while j <= ub/grid;

	print "     Grid # = " j;

	st0 = ub-j*grid;
	bt = ((alpha^2)*intk2*st0^(2*q+1)/((2*q+1)*intx2k2))^(1/(4*q+1));
	{rn} = bartlett(h,bt,q);
	{rd} = bartlett(h,bt,0);
	r2 = (rn/rd)^2;
	f0 = st0-(q*(kq^2)*r2*t/intk2)^(1/(2*q+1));

	if f0 == 0;	@ if exact solution is found @

		print "---> Grid search found an exact root!";
		print "";
		sthat = st0;
		goto next;

	elseif fub*f0 < 0;	@ if a solution is located near here @

		print "---> Grid search found a root near here!";
        print "";
		print "Bisection search starts here...";
		print "";
		stmax = st0+grid;   @ upper bound for bisection search @
		stmin = st0;        @ lower bound for bisection search @
		goto bisect;

	else;

		goto cont;

	endif;

	cont:

j = j+1; endo;

@-----------------------------------------------------------------------------@
@   SECOND, BISECTION SEARCH ROUTINE IS GIVEN BELOW.                          @
@-----------------------------------------------------------------------------@

bisect:

j = 1; do until abs(stmax-stmin) < maxtol;

	print "     Iteration # =" j;

	stmid = (stmin+stmax)/2;	@ midpoint of the interval [stmin,stmax] @

	btmax = ((alpha^2)*intk2*stmax^(2*q+1)/((2*q+1)*intx2k2))^(1/(4*q+1));
	{rn} = bartlett(h,btmax,q);
	{rd} = bartlett(h,btmax,0);
	r2 = (rn/rd)^2;
	fmax = stmax-(q*(kq^2)*r2*t/intk2)^(1/(2*q+1));

	btmid = ((alpha^2)*intk2*stmid^(2*q+1)/((2*q+1)*intx2k2))^(1/(4*q+1));
	{rn} = bartlett(h,btmid,q);
	{rd} = bartlett(h,btmid,0);
	r2 = (rn/rd)^2;
	fmid = stmid-(q*(kq^2)*r2*t/intk2)^(1/(2*q+1));

	if fmid == 0;	@ if an exact solution is found @

		sthat = stmid;
		goto next;

	elseif fmid*fmax < 0;	@ if a solution is between fmid and fmax @

		stmin = stmid;

	else;	@ if a solution is between fmin and fmid @	

		stmax = stmid;

	endif;

j = j+1; endo;

print "---> Bisection search found a root!";
print "";

sthat = stmax;

@-----------------------------------------------------------------------------@
@   HAC COVARIANCE ESTIMATE "OMEGA_HAT" IS OBTAINED HERE.                     @
@-----------------------------------------------------------------------------@

next:

{omegahat} = bartlett(g,sthat,0);



/****************************************************************************/;
/*****   3. Obtaining the graphical solution of optimal bandwidth       *****/;                  
/****************************************************************************/;

@-----------------------------------------------------------------------------@
@   YOU CAN CHECK THE OPTIMAL BANDWIDTH GRAPHICALLY.                          @
@-----------------------------------------------------------------------------@

s = seqa(0,1/100,5000);  @ ADJUST THE NUMBER OF GRIDS for THE GRAPH. @
lhs = s;    @ a 45-degree line @
rhs = s;

j = 1; do while j <= rows(rhs);

	bt = ((alpha^2)*intk2*(s[j])^(2*q+1)/((2*q+1)*intx2k2))^(1/(4*q+1));
	{rn} = bartlett(h,bt,q);
	{rd} = bartlett(h,bt,0);
	r2 = (rn/rd)^2;
	rhs[j] = (q*(kq^2)*r2*t/intk2)^(1/(2*q+1));

j = j+1; endo;

library pgraph;
graphset;
fonts("complex");
title("Graphical Solution of SP Bandwidth for Bartlett HAC Estimator");
xlabel("S(T)_hat");
_pltype = {4,6};    @ line type @
_pcolor = {1,2};    @ line color @
_plwidth = {10,10}; @ line width @
xy(s,lhs~rhs);



/****************************************************************************/;
/*****   4. Displaying the results                                      *****/;                  
/****************************************************************************/;

@-----------------------------------------------------------------------------@
@   DEFINE YOUR OUTPUT FILE HERE.                                             @
@-----------------------------------------------------------------------------@

output file = c:\iphac\gauss\bt_sp.out reset;

print "******************************************************************************";
print "******************************************************************************";
print "*******   HAC COVARIANCE ESTIMATE WITH BARTLETT KERNEL & SP BANDWIDTH   ******";
print "******************************************************************************";
print "******************************************************************************";
print "";
print "   DATE = " datestr(0);;
print "   TIME = " timestr(0);
print "";
print "   SAMPLE SIZE = " t;;
print "   BANDWIDTH = " sthat;
print "";
print "******************************************************************************";
print "******************************************************************************";
format /m1 /ldn 8,4;
print "";
print omegahat;
print "";
print "******************************************************************************";
print "******************************************************************************";



/****************************************************************************/;
/*****   A. Subroutine: "proc bartlett"                                 *****/;
/*****   Estimation of the generalized derivative of spectral density   *****/;
/*****   at the origin with the Bartlett kernel                         *****/;                  
/****************************************************************************/;

proc(1) = bartlett(g,st,q);
local maxlagbt,jbt,sbt,s0bt,xbt,wbt;

	if q == 0; sbt=(g'g)/t;
	else; sbt = zeros(cols(g),cols(g));
	endif;

	if st <= t-1; maxlagbt = st;
	else; maxlagbt = t-1;
	endif;

	jbt = 1; do while jbt <= maxlagbt;

		xbt = jbt/maxlagbt;
		s0bt=(g[1:t-jbt,.]'g[1+jbt:t,.])/t;

		if xbt <= 1; wbt = 1-xbt;
		else; wbt = 0;
		endif;

		sbt = sbt+wbt*(jbt^q)*(s0bt+s0bt');

	jbt = jbt+1; endo;

retp(sbt);
endp;