*
* Procedure for calculating Hurst exponent for a single series
*
* Source    : "Chaos and Order in the Capital Markets," by Edgar Peters.
*             John Wiley & Sons, Inc., New York, 1991
*
* Syntax:
*
*   @HURST(options)  series  start  end
*
*   Parameters:
*      series     --  the series to analyze
*      start,end  --  specifies the range to graph, can be omitted
*
*   Options:
*      GRAPH/[NOGRAPH]  -- specifies whether you want to plot graphs of
*                          LOG(R/S) vs. LOG(N).
*      HEADER=header for graph
*
*   This procedure will:
*     1) Store the calculated LOG(R/S) and LOG(# of OBS) to a PRN
*        format file called HURST.PRN
*     2) Plot the LOG(R/S) vs. LOG(# of OBS), (base 10 logs,
*        not natural logs).  Only if you use GRAPH option
*     3) Prompt you to select the range over which you want to run
*        the regression of LOG(R/S) on CONSTANT and LOG(# of OBS) to
*        get the Hurst exponent (coefficient on the LOG(# of obs.) series).
*     4) Do the regression, storing the calculated hurst exponent in
*        a global real variable called %%HURST.
*     5) Repeat the first plot, this time with a line representing
*        the calculated Hurst exponent.
*
* Revision Schedule:
*    07/93 Written by Tom Maycock and Tom Doan, Estima
*    12/06 Revised to use DBOX rather than multiple queries. Improved
*          labelling on graphs.
*
procedure hurst series start end
type series series
type integer start end
*
option switch graph  1
option string header
*
local integer nbeg nend nobs mend n m count
local real lbeg_reg lend_reg
local integer beg_reg end_reg
local series ave_range dev total numobs logrs lognobs hser range
*
declare real %%hurst
*
inquire(series=series) nbeg<<start nend<<end
compute nobs=nend-nbeg+1
clear range ave_range dev
smpl nbeg nend
*
*  Loop over range of values for N
*
do n=1,nobs
   compute count = 1
   *
   *  Determine how many sets of N observations are in the data set.
   *  Loop over that number of blocks of data
   *
   compute mend = (((nend-nbeg+1)/n)*n)+(nbeg-1)-(n-1)
   do m=nbeg,mend,n
      *
      *  Look at current block of N observations and compute
      *  R/S range for this block
      *
      accum series m m+n-1 total
      set   dev m m+n-1 = series{0} - (total(m+n-1)/n)
      accum dev m m+n-1 dev
      extremum(noprint) dev m m+n-1
      stats(noprint) series m m+n-1
      set range count count =  (%maximum - %minimum)/sqrt(%variance)
      compute count = count + 1
   end do m
   *
   *  Now take average R/S of the various blocks
   *
   accum range 1 (count-1)
   set ave_range n n = range(count-1)/(count-1)
end do n
*
*  Generate LOG(base 10) of NOBS, and take LOG(base 10) of R/S
*
set numobs  1 nobs = t
set logrs   1 nobs = log(ave_range)/log(10.)
set lognobs 1 nobs = log(numobs)/log(10.)
*
*  Write out the results
*
*open copy hurst.prn
*copy(for=prn,org=obs) 1 nobs logrs lognobs
*
*  Plot Log(R/S) versus Log(# of obs.)
*
if (graph)
{
   scatter(vlabel="Log(R/S)",hlabel="LOG(# of observations)", $
        header=header,style=line) 1
   # lognobs logrs 1 nobs
}
*
*  Get entry range for Hurst exponent calculation
*
compute [real] lbeg_reg = lognobs(nobs)*.25
compute [real] lend_reg = lognobs(nobs)*.75

dbox(action=define,title="Regression Range")
dbox(stext="Enter lower and upper bound for regressions in")
dbox(stext="terms of log(# of observations)")
dbox(stext="Lower",width=10)
dbox(edittext,special=sameline,init,width=12) lbeg_reg
dbox(stext="Upper",width=10)
dbox(edittext,special=sameline,init,width=12) lend_reg
dbox(action=run)

compute beg_reg = fix(exp(lbeg_reg*log(10.)))
compute end_reg = fix(exp(lend_reg*log(10.)))
*
*  Calculate Hurst Exponent
*
linreg logrs beg_reg end_reg
# constant lognobs
compute %%hurst = %beta(2)
display "Hurst exponent = " %%hurst
if (graph)
{
   *
   *  Calculate fitted Hurst exponent, and plot
   *
   prj hser beg_reg end_reg
   *
   spgraph
   scatter(vlabel="Log(R/S)",hlabel="LOG(# of observations)",$
     style=line,header=header,lines=||%beta(1),0.5|%beta(1),1.0||) 2
   # lognobs logrs 1 nobs
   # lognobs hser  beg_reg end_reg
   grtext(position=upleft) "Hurst exponent = "+%strval(%%hurst,"*.##")
   spgraph(done)
}
end hurst









