*! whittle Stan Hurn/Ken Lindsay/Kit Baum *! performs local Whittle and exact local Whittle estimation of the memory parameter d *! 1.0.1 15apr2019 *! 1.0.2 26apr2019 remove adjust option *! 1.0.3 20jun2019 add logic from gphudak *! 1.0.4 01nov2019 ensure that d=1 is not passed to fracdiff *! 1.0.5 06mar2023 return scalars to enable rolling capt prog drop whittle mata:mata clear prog whittle, rclass byable(recall) version 13 syntax varlist(numeric ts max=1) [if] [in] [, Powers(numlist >0 <1) DETRend EXACT] loc qq qui marksample touse // get time variables _ts timevar panelvar if `touse', sort onepanel markout `touse' `timevar' qui tsreport if `touse' if r(N_gaps) { di as err "sample may not contain gaps" exit 498 } loc st = string(`r(start)', "`r(tsfmt)'") loc nd = string(`r(end)', "`r(tsfmt)'") quietly count if `touse' if `r(N)' == 0 error 2000 loc en = `r(N)' // if powers not specified set default as 0.65 if "`powers'" == "" { loc powers = 0.65 } // do the appropriate detrending tempvar xdt if "`detrend'" == "detrend" { `qq' reg `varlist' `timevar' if `touse' loc dt "`varlist', detrended" } else { `qq' reg `varlist' if `touse' loc dt "`varlist'" } qui predict double `xdt' if e(sample), res // set up the format of the reporting loc Npowers: word count `powers' tempname whit mat `whit' = J(`Npowers', 6, 0) mat colnames `whit' = N power trunc d se ase di as res _dup(60) "-" di as res " N Power Trunc Est d StdErr Asy.StdErr" di as res _dup(60) "-" // call the routines tokenize `powers' loc i 1 while "``i''" != "" { loc trunc = floor(`en'^``i'') loc eye = ``i'' if "`exact'" == "" { mata: calcwhittle("`xdt'","`touse'",`trunc',"`st'","`nd'","`dt'") } else { mata: calcewhittle("`xdt'","`touse'",`trunc',"`st'","`nd'","`dt'") } di as res %6.0f `en' " " %5.2f `eye' " " %5.0f `trunc' /// " " %8.0g d__ " " %8.4g se__ " " %8.4g ase__ " mat `whit'[`i',1] = `en', ``i'', `trunc', d__, se__, ase__ loc i = `i' + 1 } di as res _dup(60) "-" return local cmdname "whittle" return local varname `varlist' mat `whit' = `whit'' ret scalar N = `whit'[1,1] ret scalar power = `whit'[2,1] ret scalar trunc = `whit'[3,1] ret scalar d = `whit'[4,1] ret scalar se = `whit'[5,1] ret scalar ase = `whit'[6,1] ret matrix whittle `whit' end version 13 mata: void function calcwhittle( string scalar xdt, string scalar touse, real scalar m, string scalar st, string scalar nd, string scalar dt) { st_view(x=.,.,xdt,touse) n=rows(x) // call routine to do local Whittle estimation el = whittle(x,m) // pass results to Stata st_numscalar("N",n) st_numscalar("d__",el[1,1]) st_numscalar("se__",el[1,2]) st_numscalar("ase__",sqrt(0.25)/sqrt(m)) } real vector whittle(x,m) { /* --------------------------------------------------------------------------- Returns the estimated value of the fractional difference parameter, d. INPUT x: data (n*1 vector) m: truncation number OUTPUT d: estimated fractional difference parameter This program only works for a single time series. Estimation is by Golden Section Search. Search is hardwired to the interval [-0.5,1.0] ----------------------------------------------------------------------------- */ n = rows(x) // computing the Fourier coefficients tn = 1..n tm = 1..m fac = 2*pi()/n arg = J(m,n,fac) tt = tm'*tn arg = arg :* tt tc = cos(arg) ts = sin(arg) wxr = (tc * x) :/ n wxi = -(ts * x) :/ n // power spectrum Ix = (wxr:^2 + wxi:^2) // Fourier angular frequencies lambda = fac :* tn' Ix = Ix[1..m] lambda = lambda[1..m] // Golden Section Search - pre-assigned values of xl,xr fixes search interval to (-0.5,1.0) xl = -0.5 xr = 1.0 gratio = 0.618033988749890 tol = 5.0e-7 llm = mean(log(lambda)) xlower = xl+(xr-xl)*gratio^2 vlower = log( mean(lambda:^(2*xlower) :* Ix )) - 2*xlower*llm xupper = xl+(xr-xl)*gratio vupper = log( mean(lambda:^(2*xupper) :* Ix )) - 2*xupper*llm do { if ( vlower>=vupper ) { xl = xlower xlower = xupper vlower = vupper xupper = xl+(xr-xl)*gratio vupper = log( mean(lambda:^(2*xupper) :* Ix )) - 2*xupper*llm } else { xr = xupper; xupper = xlower vupper = vlower xlower = xl+(xr-xl)*gratio^2 vlower = log( mean(lambda:^(2*xlower) :* Ix )) - 2*xlower*llm } } while ( xr-xl>tol ) dnew = 0.5*(xr+xl) // test boundaries flag = 1 if ( (dnew+0.5) < 5.0*tol ) { printf("\nUnreliable Whittle estimate - value too close to lower boundary\n") flag = 0 } if ( (1.0-dnew) < 5.0*tol ) { printf("\nUnreliable Whittle estimate - value too close to upper boundary\n") flag = 0 } // set up values to return if (flag) { d0 = mean( lambda:^(2*dnew) :* Ix ) d1 = 2*mean( log(lambda) :* lambda:^(2*dnew) :* Ix ) d2 = 4*mean( (log(lambda)):^2 :* lambda:^(2*dnew) :* Ix ) se = d0/(sqrt(m)*sqrt(d0*d2-d1^2)) x2 = (dnew,se ) return(x2 ) } else { return( dnew , -999 ) } } // --------------------------------------------------------------------------- void function calcewhittle( string scalar xdt, string scalar touse, real scalar m, string scalar st, string scalar nd, string scalar dt) { st_view(x=.,.,xdt,touse) n=rows(x) // call routine to do exact local Whittle estimation el3 = ewhittle(x,m) // pass results back to Stata st_numscalar("N",n) st_numscalar("d__",el3[1,1]) st_numscalar("se__",el3[1,2]) st_numscalar("ase__",sqrt(0.25)/sqrt(m)) } /*----------------------------------------------------------------------------- Compute the exact Whittle likelihood Shimotsu and Phillips (2004) and Shimotsu (2004). INPUT x: data (n*1 vector) d: fractional index m: truncation number -----------------------------------------------------------------------------*/ real scalar function ewfeval(x,d,m) { n=rows(x) dx=fracdiff(x,d) // calling built-in Mata function // computing the Fourier coefficients tn = 1..n tm = 1..m fac = 2*pi()/ n arg = J(m,n,fac) tt = tm'*tn arg = arg :* tt tc = cos(arg) ts = sin(arg) wxr = (tc * dx) :/ n wxi = -(ts * dx) :/ n // power spectrum Ix = (wxr:^2 + wxi:^2) // Fourier angular frequencies lambda = fac :* tn' Ix = Ix[1..m] lambda = lambda[1..m] // return function value for exact local Whittle likelihood g = mean(Ix) r = log(g) - 2*d*mean(log(lambda)) return(r) } /*----------------------------------------------------------------------------- Compute the fractional difference index based on the exact Whittle likelihood of Shimotsu and Phillips (2004) and Shimotsu (2004). INPUT x: data (n*1 vector) m: truncation number -----------------------------------------------------------------------------*/ real vector ewhittle(x,m) { // Golden Section Search - pre-assigned values of xl,xr fixes search interval to (-0.5,1.0) xl = -0.5 xr = 1.0 gratio = 0.618033988749890 tol = 5.0e-7 xlower = xl+(xr-xl)*gratio^2 vlower = ewfeval(x,xlower,m) xupper = xl+(xr-xl)*gratio vupper = ewfeval(x,xupper,m) do { if ( vlower>=vupper ) { xl = xlower xlower = xupper vlower = vupper xupper = xl+(xr-xl)*gratio vupper = ewfeval(x,xupper,m) } else { xr = xupper xupper = xlower vupper = vlower xlower = xl+(xr-xl)*gratio^2 vlower = ewfeval(x,xlower,m) } } while ( xr-xl>tol ) dnew = 0.5*(xr+xl) // test boundaries flag = 1 if ( (dnew+0.5) < 0.01 ) { printf("\nUnreliable eWhittle estimate - value too close to lower boundary\n") flag = 0 } if ( (1.0-dnew) < 0.01 ) { printf("\nUnreliable eWhittle estimate - value too close to upper boundary\n") flag = 0 } // finite difference estimate of second derivative if (flag) { dl = dnew*0.99 du = dnew*1.01 fl = ewfeval(x,dl,m) fc = ewfeval(x,dnew,m) fu = ewfeval(x,du,m) d2 = 1.0e4*(fl - 2*fc + fu)/dnew^2 // Test for convexity and return if ( d2 > 0 ) { x2 = ( dnew, sqrt(1/(m*d2)) ) return( x2 ) } else { printf("\neWhittle likelihood function is not convex\n"); return( dnew , 0 ) } } else { return( dnew , -999 ) } } end