*
* @HURST(options)  series  start  end
* Calculates the Hurst exponent for a single series
*
* Parameters:
*   series     --  the series to analyze
*   start,end  --  specifies the range to graph, can be omitted
*
* Options:
*   [GRAPH]/NOGRAPH
*     Controls graphing of LOG(R/S) vs. LOG(N).
*   HEADER=header for graph
*
*   [PRINT]/NOPRINT
*      Controls printing of the log-log regression
*
*   [DIALOG]/NODIALOG
*     Controls whether you are prompted for the range over which to
*     run the log-log regression
*
* This procedure will:
*   1) Plot the LOG(R/S) vs. LOG(# of OBS), (base 10 logs,
*      not natural logs).  Only if you use GRAPH option
*   2) 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).
*      (If you use NODIALOG, it will use a standard range)
*   3) Do the regression, storing the calculated hurst exponent in
*      a global real variable called %%HURST.
*   4) Repeat the first plot, this time with a line representing
*      the calculated Hurst exponent.
*
* Variables Defined:
*   %%HURST   Hurst exponent
*
* Revision Schedule:
*   07/1993 Written by Tom Maycock and Tom Doan, Estima
*   12/2006 Revised to use DBOX rather than multiple queries. Improved
*           labelling on graphs.
*   08/2010 Add infobox
*   03/2014 Added DIALOG/NODIALOG and PRINT/NOPRINT options
*   03/2014 Removed unnecessary SMPL instruction
*
procedure hurst series start end
type series series
type integer start end
*
option switch graph  1
option string header
option switch dialog 1
option switch print  1
*
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
set ave_range 1 nobs = 0
set range     1 nobs = 0
*
*  Loop over range of values for N
*
infobox(action=define,progress,lower=1,upper=nobs) "Calculating R/S Statistics"
do n=3,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
      *
      diff(center) series m m+n-1 dev
      accum dev m m+n-1 dev
      stats(noprint) series m m+n-1
      extremum(noprint) dev m m+n-1
      compute range(count)=(%maximum - %minimum)/sqrt(%variance)
      compute count = count + 1
   end do m
   *
   *  Now take average R/S of the various blocks
   *
   sstats(mean) 1 count-1 range>>ave_range(n)
   infobox(current=n)
end do n
infobox(action=remove)
*
*  Generate LOG(base 10) of NOBS, and take LOG(base 10) of R/S
*
set numobs  3 nobs = t
set logrs   3 nobs = log(ave_range)/log(10.)
set lognobs 3 nobs = log(numobs)/log(10.)
*
*  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 3 nobs
}
*
*  Get entry range for Hurst exponent calculation
*
compute [real] lbeg_reg = lognobs(nobs)*.25
compute [real] lend_reg = lognobs(nobs)*.75

if dialog {
   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(print=print,title="Estimation of Hurst Exponent") logrs beg_reg end_reg
# constant lognobs
compute %%hurst = %beta(2)
if print
   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 3 nobs
   # lognobs hser  beg_reg end_reg
   grtext(position=upleft) "Hurst exponent = "+%strval(%%hurst,"*.##")
   spgraph(done)
}
end hurst
