[R-SIG-Finance] TTR Yang-Zhang volatility bug?

rex rex at nosyntax.net
Fri Mar 15 20:49:14 CET 2013


I'm generating lognormal OHLC series and testing the RNG using TTR's
volatility functions. The yang.zhang calculation prepends too many NAs.
When Joshua's Yang-Zhang code is added as a function in my code, it
prepends fewer NAs, but still one too many, I think. 

Why does what appears to be the same code behave differently?

R version 2.15.3 (2013-03-01) -- "Security Blanket"
Copyright (C) 2013 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-pc-linux-gnu (64-bit)

   library(TTR)
   library(xts)

   vYZ = function(OHLC, n=10, N=252, ...){     #Joshua's Yang-Zhang TTR function
   # Historical Open-High-Low-Close Volatility: Yang Zhang
      # http://www.sitmo.com/eq/417            #dead link
      if(is.xts(OHLC)) {
        Cl1 <- lag(OHLC[,4])
      } else {
        Cl1 <- c( NA, OHLC[-NROW(OHLC),4] )
      }
      
      dots <- list(...)
      if(is.null(dots$alpha)) {
        alpha <- 1.34
      }
      if(is.null(dots$k)) {
        k <- (alpha-1) / ( alpha + (n+1)/(n-1) )
      }

      s2o <- N * runVar(log(OHLC[,1] / Cl1), n=n)
      s2c <- N * runVar(log(OHLC[,4] / OHLC[,1]), n=n)
      #s2rs <- volatility(OHLC=OHLC, n=n, calc="rogers.satchell", N=N, ...)
      s2rs <- sqrt( N/n * runSum(        #taken from Joshua's rogers.satchell code
                 log(OHLC[,2]/OHLC[,4]) * log(OHLC[,2]/OHLC[,1]) +
                 log(OHLC[,3]/OHLC[,4]) * log(OHLC[,3]/OHLC[,1]), n ) )

      s <- sqrt(s2o + k*s2c + (1-k)*(s2rs^2))
   }

   #my code starts here
   options(width=132, digits=4)
   ohlc           = as.zoo(matrix(nrow=nDays,ncol=4))
   colnames(ohlc) = c('Open','High','Low','Close')
   allDays        = seq(as.Date('2012-01-01'), as.Date('2012-02-24'), by = 1)
   wkDays         = allDays[!weekdays(allDays) %in% c('Saturday','Sunday')]
   index(ohlc)    = wkDays
   ohlc           = as.xts(ohlc)
   tDays          = 252                         #trading days/year
   nDays          = length(wkDays)              #random walk this many trading days
   hist           = 15                          #days used to estimate volatility
   drift          = 0.10                        #annual trend of walk
   vol            = 0.25                        #25% annual volatility
   stkPr          = 100                         #initial stock price
   pricesDay      = 6.5*12                      #6.5 hours open and 12 prices per hour
   pricesWlk      = pricesDay*nDays             #number of prices in each random walk
   pricesYr       = pricesDay*tDays             #number of prices in a year
   nWalks         = 100                         #number of random walks
   tRands         = pricesWlk*nWalks            #number of lognormal rands required
   mu             = drift/(pricesYr)            #trend per trade for walk
   sd             = vol/sqrt(pricesYr)          #volatility per trade for walk          
   expPr          = stkPr*exp((drift + 0.5*vol*vol)*nDays/tDays)      #expected end price
   y              = matrix(rlnorm(tRands, mu, sd), nWalks, pricesWlk) #random lognormals
   z              = stkPr*t(apply(y,1,cumprod)) #do the day-by-day multiplication
   endPr          = mean(z[,pricesWlk])         #average price at end of walk
   
   sumYZ = sumCC = 0
   for (w in 1:nWalks){
     for (i in 1:nDays){
       j         = (i-1)*pricesDay + 1
       k         = j + pricesDay -1
       ohlc[i,1] = z[w,j]           #open for day
       ohlc[i,2] = max(z[w,j:k])    #high for day
       ohlc[i,3] = min(z[w,j:k])    #low for day
       ohlc[i,4] = z[w,i*pricesDay] #close for day
     }
     volYZ2 = vYZ(ohlc, n=hist, N=tDays)
     volYZ = volatility(ohlc, n=hist, N=tDays, calc="yang.zhang") #Yang-Zhang volatility
     volCC = volatility(ohlc, n=hist, N=tDays, calc="close")      #closing volatility
     sumYZ = sumYZ + mean(na.omit(volYZ))
     sumCC = sumCC + mean(na.omit(volCC))
     #print(c(mean(volYZ),mean(volCC)))     #there are multiple estimates
   }
   print(c('averageYZ', sumYZ/nWalks))
   #[1] "averageYZ"         "0.227380100996657"
   print(c('averageCC', sumCC/nWalks))
   #[1] "averageCC"         "0.239352679348972"
   #chartSeries(ohlc, name='SIM', theme='white')

> volYZ             #too many NAs
               [,1]
2012-01-02     NA
2012-01-03     NA
2012-01-04     NA
2012-01-05     NA
2012-01-06     NA
2012-01-09     NA
2012-01-10     NA
2012-01-11     NA
2012-01-12     NA
2012-01-13     NA
2012-01-16     NA
2012-01-17     NA
2012-01-18     NA
2012-01-19     NA
2012-01-20     NA
2012-01-23     NA
2012-01-24     NA
2012-01-25     NA
2012-01-26     NA
2012-01-27     NA
2012-01-30     NA
2012-01-31     NA
2012-02-01     NA
2012-02-02     NA
2012-02-03     NA
2012-02-06     NA
2012-02-07     NA
2012-02-08     NA
2012-02-09     NA
2012-02-10 0.2125
2012-02-13 0.2042
2012-02-14 0.2010
2012-02-15 0.2011
2012-02-16 0.2077
2012-02-17 0.2133
2012-02-20 0.2222
2012-02-21 0.2261
2012-02-22 0.2270
2012-02-23 0.2233
2012-02-24 0.2271

> volYZ2             #better
               [,1]
2012-01-02     NA
2012-01-03     NA
2012-01-04     NA
2012-01-05     NA
2012-01-06     NA
2012-01-09     NA
2012-01-10     NA
2012-01-11     NA
2012-01-12     NA
2012-01-13     NA
2012-01-16     NA
2012-01-17     NA
2012-01-18     NA
2012-01-19     NA
2012-01-20     NA
2012-01-23 0.2067
2012-01-24 0.2087
2012-01-25 0.2100
2012-01-26 0.2110
2012-01-27 0.2139
2012-01-30 0.2138
2012-01-31 0.2091
2012-02-01 0.2163
2012-02-02 0.2240
2012-02-03 0.2272
2012-02-06 0.2245
2012-02-07 0.2259
2012-02-08 0.2229
2012-02-09 0.2183
2012-02-10 0.2151
2012-02-13 0.2069
2012-02-14 0.2038
2012-02-15 0.2038
2012-02-16 0.2099
2012-02-17 0.2155
2012-02-20 0.2229
2012-02-21 0.2266
2012-02-22 0.2292
2012-02-23 0.2250
2012-02-24 0.2285

> volCC
               [,1]
2012-01-02     NA
2012-01-03     NA
2012-01-04     NA
2012-01-05     NA
2012-01-06     NA
2012-01-09     NA
2012-01-10     NA
2012-01-11     NA
2012-01-12     NA
2012-01-13     NA
2012-01-16     NA
2012-01-17     NA
2012-01-18     NA
2012-01-19     NA
2012-01-20 0.2288
2012-01-23 0.2281
2012-01-24 0.2628
2012-01-25 0.2618
2012-01-26 0.2715
2012-01-27 0.2710
2012-01-30 0.2766
2012-01-31 0.2603
2012-02-01 0.2891
2012-02-02 0.3014
2012-02-03 0.3083
2012-02-06 0.2822
2012-02-07 0.2812
2012-02-08 0.2813
2012-02-09 0.2788
2012-02-10 0.2845
2012-02-13 0.2471
2012-02-14 0.2654
2012-02-15 0.2560
2012-02-16 0.2827
2012-02-17 0.2627
2012-02-20 0.2832
2012-02-21 0.2686
2012-02-22 0.3048
2012-02-23 0.3093
2012-02-24 0.3536

Thanks,

-rex
--



More information about the R-SIG-Finance mailing list