/***********************************************************************************/
/*   Program to apply modified GPH estimator for long memory estimation            */
/*   See "Level Shifts and the Illusion of Long Memory in Economic Time Series,"   */
/*   by Aaron Smith, JBES, 2005, 23(3):321-335.                                    */
/***********************************************************************************/

/***********************************************************************************/
/*  Note that this program will give slightly different modified-GPH estimates     */
/*  for these data than in the published paper. For k=3 in Table 4, the paper has  */
/*  0.48 and 0.42 as the estimates of the long memory parameter for the plug-in    */
/*  and root-T methods of choosing J, respectively. This program will give the     */
/*  values 0.42 and 0.55 for the corresponding cases. I made a small programming   */
/*  error in an earlier version of the code.                                       */
/***********************************************************************************/


/*  If you use this code on your own data, just give it the variable y (e.g., see  */ 
/*  line 24 below) and the values of J and k (see lines 43 and 44 below).          */


/* Read in data and create the variable y     */
loadm sandp[]=c:\research\published\longm\sandp.dat;
sp = reshape(sandp,rows(sandp)/3,3);
y=abs(250*ln(sp[2:rows(sp),2]./sp[1:rows(sp)-1,2])); 


  /* Generate plug-in value of J */
  samp=rows(y);
  t=seqa(0,1,samp);
  LL=0.1*samp^(6/7);
  lam=2*pi*seqa(1,1,LL)/samp;
  reg=ones(LL,1)~-ln(2-2*cos(lam))~(-0.5*lam^2);
  lper = ln((sumc(y.*cos(t*lam'))^2 + sumc(y.*sin(t*lam'))^2)/(2*pi*samp));
  bhat = inv(reg'reg)*reg'lper; 
  plug_in=floor(((27/(128*pi^2))^(1/5))*(maxc(bhat[3]|1)^(-2/5))*(samp^(4/5))); 


/***************************************************************/
/*   Choose J and k                                            */
/*     For plug_in selection of J, the value of k (denoted kc) */
/*     must be an integer between 1 and 5                      */
/***************************************************************/
  J=plug_in|floor(sqrt(samp)); 
  kc=3; 
/**************************************************************/


"--------------------------------------------------------------";
format /rz 8,6;
"k = " kc; 
"T = " samp;  
format /rz 10,4;

jc = 1;
do while jc le rows(J); 
  gtt = J[jc];

  /* Generate log periodogram (lper), regressor for GPH regression (XX), and Fourier frequencies (lam) */
  lam=2*pi*seqa(1,1,gtt)/samp;
  lper = ln((sumc(y.*cos(t*lam'))^2 + sumc(y.*sin(t*lam'))^2)/(2*pi*samp));
  XX=ones(gtt,1)~-ln(2-2*cos(lam));

  /* GPH regression  */
  bhat = inv(XX'XX)*XX'lper;
  cse = sqrt(diag(inv(XX'XX)*((pi^2)/6)));

  /* Generate appropriate J for modified GPH regression (plug-in method requires that J be scaled up depending on k)  */
  if jc eq 1;
    if kc eq 1;      gtt1=floor(1.69*gtt);
    elseif kc eq 2;  gtt1=floor(1.88*gtt);
    elseif kc eq 3;  gtt1=floor(2.09*gtt);
    elseif kc eq 4;  gtt1=floor(2.33*gtt);
    elseif kc eq 5;  gtt1=floor(2.58*gtt);
    endif;
    ""; "                  Plug-in selection of J";
  else;
    gtt1=gtt;
    ""; "                         Fixed J";
  endif;

  /* Generate log periodogram (lper), regressor for mod-GPH regression (zz), and Fourier frequencies (lam) */
  lam=2*pi*seqa(1,1,gtt1)/samp;
  lper = ln((sumc(y.*cos(t*lam'))^2 + sumc(y.*sin(t*lam'))^2)/(2*pi*samp));
  k=kc*gtt/samp;
  zz = ones(gtt1,1)~-ln(2-2*cos(lam))~-ln(k^2+lam^2);

  /* modified GPH regression  */
  cbhat = inv(zz'zz)*zz'lper;
  ccse = sqrt(diag(inv(zz'zz)*((pi^2)/6)));

/* Print results */
"                 J         d        s.e.       t-stat";
"GPH:     " gtt~bhat[2]~cse[2]~(bhat[2]/cse[2]);
"mod GPH: " gtt1~cbhat[2]~ccse[2]~(cbhat[2]/ccse[2]); 

jc=jc+1;
endo;
"--------------------------------------------------------------";