/****************************************************************************/;
/****************************************************************************/;
/*****   Solve-the-Equation Plug-In (SP) Bandwidth Choice Rule for      *****/;
/*****   HAC Estimation with the Parzen 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 = 2;                  @ characteristic exponent of Parzen kernel @			
kq = 6;                 @ generalized derivative of Parzen kernel at origin @
intk2 = 151/280;        @ squared integral of Parzen kernel @
intx4k2 = 929/295680;   @ "4th-order moment" of Bartlett kernel @



/****************************************************************************/;
/*****   2. Implementing Parzen-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 = -(1+8*phi+phi^2)/(phi^2-1)^2;	@ 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)*intx4k2))^(1/(4*q+1));
{rn} = parzen(h,bt,q);    @ estimate of q'th generalized derivative of spectral density at origin @
{rd} = parzen(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)*intx4k2))^(1/(4*q+1));
	{rn} = parzen(h,bt,q);
	{rd} = parzen(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)*intx4k2))^(1/(4*q+1));
	{rn} = parzen(h,btmax,q);
	{rd} = parzen(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)*intx4k2))^(1/(4*q+1));
	{rn} = parzen(h,btmid,q);
	{rd} = parzen(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} = parzen(g,sthat,0);



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

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

s = seqa(0,1/100,2000);  @ 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)*intx4k2))^(1/(4*q+1));
	{rn} = parzen(h,bt,q);
	{rd} = parzen(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 Parzen 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\pz_sp.out reset;

print "******************************************************************************";
print "******************************************************************************";
print "********   HAC COVARIANCE ESTIMATE WITH PARZEN 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 parzen"                                   *****/;
/*****   Estimation of the generalized derivative of spectral density   *****/;
/*****   at the origin with the Parzen kernel                           *****/;                  
/****************************************************************************/;

proc(1) = parzen(g,st,q);
local  maxlagpz,jpz,spz,s0pz,xpz,wpz;

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

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

	jpz = 1; do while jpz <= maxlagpz;

		xpz = jpz/maxlagpz;
		s0pz = (g[1:t-jpz,.]'g[1+jpz:t,.])/t;

		if xpz <= 1/2; wpz = 1-6*xpz^2+6*xpz^3;
		elseif xpz <= 1; wpz = 2*(1-xpz)^3;
		else; wpz = 0;
		endif;

		spz = spz+wpz*(jpz^q)*(s0pz+s0pz');

	jpz = jpz+1; endo;

retp(spz);
endp;