[Blotter-commits] r574 - in pkg/RTAQ: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Mar 13 19:45:47 CET 2011
Author: jonathan
Date: 2011-03-13 19:45:47 +0100 (Sun, 13 Mar 2011)
New Revision: 574
Modified:
pkg/RTAQ/R/periodicityTAQ.R
pkg/RTAQ/man/spotVol.rd
Log:
spotvol update
Modified: pkg/RTAQ/R/periodicityTAQ.R
===================================================================
--- pkg/RTAQ/R/periodicityTAQ.R 2011-03-12 18:52:40 UTC (rev 573)
+++ pkg/RTAQ/R/periodicityTAQ.R 2011-03-13 18:45:47 UTC (rev 574)
@@ -1,51 +1,66 @@
-
-
# Documented function:
-spotVol = function (pdata, dailyvol = "bipower" , periodicvol = "TML" , on = "minutes" , k = 5 ,
- dummies = FALSE , P1 = 4 , P2 = 2, ...)
+spotVol = function(pdata, dailyvol = "bipower", periodicvol = "TML", on = "minutes",
+ k = 5, dummies = FALSE, P1 = 4, P2 = 2, marketopen = "09:30:00",
+ marketclose = "16:00:00")
{
- dates = unique(format(time(pdata),"%Y-%m-%d"))
- cDays = length(dates);
- rdata = mR = c();
- for(d in 1:cDays){
- pdatad = pdata[as.character(dates[d])]
- pdatad = aggregatePrice(pdatad, on=on , k = k )
- rdatad = makeReturns(pdatad)
- # remove the first zero returns
- rdatad = rdatad[time(rdatad)>min(time(rdatad))]
- rdata = rbind( rdata , rdatad )
- mR = rbind( mR , as.numeric(rdatad))
+ require(chron);
+ dates = unique(format(time(pdata), "%Y-%m-%d"))
+ cDays = length(dates)
+ rdata = mR = c()
+ if(on=="minutes"){
+ intraday = seq(from=times(marketopen), to=times(marketclose), by=times(paste("00:0",k,":00",sep="")))
}
- M = ncol(mR);
- if( cDays == 1){
- mR = as.numeric(rdata)
- estimdailyvol = switch( dailyvol , "bipower" = RBPCov(mR) ,
- "medrv" = MedRV(mR) ,
- "rv" = RV(mR) )
- }else{
- estimdailyvol = switch( dailyvol , "bipower" = apply(mR,1,'RBPCov') ,
- "medrv" = apply(mR,1,'MedRV') ,
- "rv" = apply(mR,1,'RV') )
+ if(tail(intraday,1)!=marketclose){intraday=c(intraday,marketclose)}
+ intraday = intraday[2:length(intraday)];
+ for (d in 1:cDays) {
+ pdatad = pdata[as.character(dates[d])]
+ pdatad = aggregatePrice(pdatad, on = on, k = k , marketopen = marketopen, marketclose = marketclose)
+ z = xts( rep(1,length(intraday)) , order.by = timeDate( paste(dates[d],as.character(intraday),sep="") , format = "%Y-%m-%d %H:%M:%S"))
+ pdatad = merge.xts( z , pdatad )$pdatad
+ pdatad = na.locf(pdatad)
+ rdatad = makeReturns(pdatad)
+ rdatad = rdatad[time(rdatad) > min(time(rdatad))]
+ rdata = rbind(rdata, rdatad)
+ mR = rbind(mR, as.numeric(rdatad))
}
- if( cDays <= 50 ){
- print("Periodicity estimation requires at least 50 observations. Periodic component set to unity")
- estimperiodicvol = rep(1,M)
- }else{
- mstdR = mR / sqrt( estimdailyvol*(1/M) );
- estimperiodicvol = diurnal( stddata = mstdR , method = periodicvol , dummies = dummies , P1 = P1 , P2 = P2 )[[1]]
- mfilteredR = mR/matrix( rep(estimperiodicvol,cDays) , byrow = T , nrow = cDays )
- estimdailyvol = switch( dailyvol , "bipower" = apply(mfilteredR,1,'RBPCov') ,
- "medrv" = apply(mfilteredR,1,'MedRV') ,
- "rv" = apply(mfilteredR,1,'RV') )
+ mR[is.na(mR)]=0
+ M = ncol(mR)
+ if (cDays == 1) {
+ mR = as.numeric(rdata)
+ estimdailyvol = switch(dailyvol, bipower = RBPCov(mR),
+ medrv = MedRV(mR), rv = RV(mR))
+ }else {
+ estimdailyvol = switch(dailyvol, bipower = apply(mR,
+ 1, "RBPCov"), medrv = apply(mR, 1, "MedRV"), rv = apply(mR,
+ 1, "RV"))
}
- out = cbind( rdata ,
- rep( sqrt( estimdailyvol*(1/M) ),each=M )*rep(estimperiodicvol,cDays) ,
- rep( sqrt( estimdailyvol*(1/M) ),each=M ) ,
- rep(estimperiodicvol,cDays) )
- out = xts( out , order.by = time(rdata) )
- names(out) = c("returns","vol","dailyvol","periodicvol")
- return( out )
+ if (cDays <= 50) {
+ print("Periodicity estimation requires at least 50 observations. Periodic component set to unity")
+ estimperiodicvol = rep(1, M)
+ }
+ else {
+ mstdR = mR/sqrt(estimdailyvol * (1/M))
+ selection = c(1:M)[ (nrow(mR)-apply(mR,2,'countzeroes')) >=20]
+ # preferably no na is between
+ selection = c( min(selection) : max(selection) )
+ mstdR = mstdR[,selection]
+ estimperiodicvol_temp = diurnal(stddata = mstdR, method = periodicvol,
+ dummies = dummies, P1 = P1, P2 = P2)[[1]]
+ estimperiodicvol = rep(1,M)
+ estimperiodicvol[selection] = estimperiodicvol_temp
+ mfilteredR = mR/matrix(rep(estimperiodicvol, cDays),
+ byrow = T, nrow = cDays)
+ estimdailyvol = switch(dailyvol, bipower = apply(mfilteredR,
+ 1, "RBPCov"), medrv = apply(mfilteredR, 1, "MedRV"),
+ rv = apply(mfilteredR, 1, "RV"))
+ }
+ out = cbind(rdata, rep(sqrt(estimdailyvol * (1/M)), each = M) *
+ rep(estimperiodicvol, cDays), rep(sqrt(estimdailyvol *
+ (1/M)), each = M), rep(estimperiodicvol, cDays))
+ out = xts(out, order.by = time(rdata))
+ names(out) = c("returns", "vol", "dailyvol", "periodicvol")
+ return(out)
}
@@ -67,115 +82,144 @@
return( 0.7413*min( M[2,]-M[1,] ) );
}
-diurnal = function( stddata , method="TML", dummies = F, P1 = 6 , P2 = 4)
+diurnal =
+function (stddata, method = "TML", dummies = F, P1 = 6, P2 = 4)
{
- # Function that has as input the standardized returns and computes the seasonality factor
- # by the scale-based or regression approach as explained in
- # Boudt, Croux and Laurent (2008)'s paper ``Robust estimation of intraweek periodicity and jump detection''
-
- # stddata is matrix that holds in row the standardized intraday returns
-
- cDays = dim(stddata)[1] ; intraT = dim(stddata)[2];
-
- meannozero = function( series ){ return(mean(series[series!=0]) ) }
- shorthscalenozero = function( series ){ return(shorthscale(series[series!=0]) ) }
- WSDnozero = function( weights , series ){
- out = sum( (weights*series^2)[series!=0] )/sum( weights[series!=0]) ;
- return( sqrt( 1.081*out) )
- }
- if(method=="SD"|method=="OLS"){ seas = sqrt(apply( stddata^2 , 2 , 'meannozero' ) ) }
- if(method=="WSD"|method=="TML"){
- seas = apply( stddata , 2 , 'shorthscalenozero' ) ;
- shorthseas = seas/sqrt(mean( seas^2) )
- shorthseas[shorthseas==0] = 1;
- weights = matrix( HRweight( as.vector(t(stddata^2) / rep(shorthseas,cDays)^2) ,qchisq(0.99,df=1)) ,
- ncol=dim(stddata)[2] , byrow=T)
- for( c in 1:intraT ){ seas[c] = WSDnozero( weights[,c] , stddata[,c] ) }
- }
- seas = seas/sqrt(mean( seas^2) )
- if( method=="OLS"|method=="TML"){
-
-
- # Do periodicity estimation on the basis of a FFF specification + Almond Polynomials
-
- #--Define variables--
- c=center();
- vstddata = as.vector(stddata); nobs = length(vstddata) ;
- vi = rep( c(1:intraT) , each = cDays ) # c(1,...,1,2,....2,3...., ,intraT)
- if(method=="TML"){ firststepresids = log(abs(vstddata))-c-log(rep(seas,each=cDays)) };
-
- # Regressors
- X = c();
- if(!dummies){
- if ( P1 > 0 ){ for( j in 1:P1 ){ X = cbind( X , cos(2*pi*j*vi/intraT) ) } } ;
- M1 = (intraT+1)/2 ; M2 = (2*intraT^2 + 3*intraT + 1)/6;
- ADD = (vi/M1 ) ; X = cbind(X,ADD);
- ADD = (vi^2/M2); X = cbind(X,ADD);
- if ( P2 > 0 ){ ADD= c(); for( j in 1:P2 ){ ADD = cbind( ADD , sin(2*pi*j*vi/intraT) ) }}; X = cbind( X , ADD ) ;
-
- #openingeffect
- opening = vi-0 ; stdopening = (vi-0)/80 ;
- almond1_opening = ( 1 - (stdopening)^3 );
- almond2_opening = ( 1 - (stdopening)^2 )*( opening);
- almond3_opening = ( 1 - (stdopening) )*( opening^2);
- X = cbind( X, almond1_opening , almond2_opening , almond3_opening ) ;
-
- #closing effect
- closing = max(vi)-vi ; stdclosing = (max(vi)-vi)/max(vi) ;
- almond1_closing = ( 1 - (stdclosing)^3 );
- almond2_closing = ( 1 - (stdclosing)^2 )*( closing);
- almond3_closing = ( 1 - (stdclosing) )*( closing^2);
- X = cbind( X, almond1_closing , almond2_closing , almond3_closing ) ;
-
- }else{
- for( d in 1:intraT){
- dummy = rep(0,intraT); dummy[d]=1; dummy = rep(dummy,each=cDays)
- X = cbind(X,dummy);
- }
- }
-
- # Remove the inliers
- selection = c(1:nobs)[vstddata!=0];
- vstddata = vstddata[selection] ; X = X[selection,]
- if(method == "TML" ){ firststepresids = firststepresids[selection]};
-
- # Dependent variable
- vy = matrix ( log( abs( vstddata )) , ncol=1 )-c;
-
- if( method == "OLS" ){
- Z = try(solve(t(X)%*%X),silent=T)
- if( inherits(Z, "try-error") ){
- print("X'X is not invertible. Switch to TML")
- }else{
- theta = solve( t(X)%*%X )%*%t(X)%*%vy ; rm(X);rm(vy);
- }
- }
- if( method == "TML" ){
- inittheta = rep(0,dim(X)[2]);
- l = -2.272; u = 1.6675;
- nonoutliers = c(1:length(vy))[( firststepresids > l) & (firststepresids <u ) ]
- truncvy = vy[nonoutliers ] ; rm(vy); truncX = X[nonoutliers,] ; rm(X);
- # MLE on the basis of the truncated dataset
- negtruncLLH = function(theta){
- res = truncvy - truncX%*%matrix(theta,ncol=1)
- return( mean( -res - c +exp(2*(res+c))/2 ) )
- }
- grnegtruncLLH = function(theta){
- res = truncvy - truncX%*%matrix(theta,ncol=1) ; dres = -truncX;
- return( apply( -dres +as.vector(exp(2*(res+c)))*dres , 2 , 'mean' ) )
- }
- est = optim(par=inittheta, fn= negtruncLLH, gr = grnegtruncLLH , method = "BFGS" ) ;
- theta = est$par;
- rm(truncX); rm(truncvy);
- }
- #--Output--
-
-
- plot( seas , main = "Non-parametric (dashed line) and parametric (full line) periodicity" , xlab="intraday period" , type="l",lty=3)
- seas = diurnalfit(theta=theta, P1=P1 , P2=P2, intraT = intraT , dummies = dummies)
- lines( seas , lty = 1 )
- return( list ( seas , theta ) )
- }else{ return( list(seas) ) }
+ cDays = dim(stddata)[1]
+ intraT = dim(stddata)[2]
+ meannozero = function(series) {
+ return(mean(series[series != 0]))
+ }
+ shorthscalenozero = function(series) {
+ return(shorthscale(series[series != 0]))
+ }
+ WSDnozero = function(weights, series) {
+ out = sum((weights * series^2)[series != 0])/sum(weights[series !=
+ 0])
+ return(sqrt(1.081 * out))
+ }
+ if (method == "SD" | method == "OLS") {
+ seas = sqrt(apply(stddata^2, 2, "meannozero"))
+ }
+ if (method == "WSD" | method == "TML") {
+ seas = apply(stddata, 2, "shorthscalenozero")
+ shorthseas = seas/sqrt(mean(seas^2))
+ shorthseas[shorthseas == 0] = 1
+ weights = matrix(HRweight(as.vector(t(stddata^2)/rep(shorthseas,
+ cDays)^2), qchisq(0.99, df = 1)), ncol = dim(stddata)[2],
+ byrow = T)
+ for (c in 1:intraT) {
+ seas[c] = WSDnozero(weights[, c], stddata[, c])
+ }
+ }
+ # new:
+ seas = na.locf(seas)
+ seas = seas/sqrt(mean(seas^2))
+ if (method == "OLS" | method == "TML") {
+ c = center()
+ vstddata = as.vector(stddata)
+ nobs = length(vstddata)
+ vi = rep(c(1:intraT), each = cDays)
+ if (method == "TML") {
+ if( length(vstddata)!= length(seas)*cDays ){ print(length(vstddata)); print(length(seas)); print(cDays)}
+ firststepresids = log(abs(vstddata)) - c - log(rep(seas,
+ each = cDays))
+ }
+ X = c()
+ if (!dummies) {
+ if (P1 > 0) {
+ for (j in 1:P1) {
+ X = cbind(X, cos(2 * pi * j * vi/intraT))
+ }
+ }
+ M1 = (intraT + 1)/2
+ M2 = (2 * intraT^2 + 3 * intraT + 1)/6
+ ADD = (vi/M1)
+ X = cbind(X, ADD)
+ ADD = (vi^2/M2)
+ X = cbind(X, ADD)
+ if (P2 > 0) {
+ ADD = c()
+ for (j in 1:P2) {
+ ADD = cbind(ADD, sin(2 * pi * j * vi/intraT))
+ }
+ }
+ X = cbind(X, ADD)
+ opening = vi - 0
+ stdopening = (vi - 0)/80
+ almond1_opening = (1 - (stdopening)^3)
+ almond2_opening = (1 - (stdopening)^2) * (opening)
+ almond3_opening = (1 - (stdopening)) * (opening^2)
+ X = cbind(X, almond1_opening, almond2_opening, almond3_opening)
+ closing = max(vi) - vi
+ stdclosing = (max(vi) - vi)/max(vi)
+ almond1_closing = (1 - (stdclosing)^3)
+ almond2_closing = (1 - (stdclosing)^2) * (closing)
+ almond3_closing = (1 - (stdclosing)) * (closing^2)
+ X = cbind(X, almond1_closing, almond2_closing, almond3_closing)
+ }
+ else {
+ for (d in 1:intraT) {
+ dummy = rep(0, intraT)
+ dummy[d] = 1
+ dummy = rep(dummy, each = cDays)
+ X = cbind(X, dummy)
+ }
+ }
+ selection = c(1:nobs)[vstddata != 0]
+ vstddata = vstddata[selection]
+ X = X[selection, ]
+ if (method == "TML") {
+ firststepresids = firststepresids[selection]
+ }
+ vy = matrix(log(abs(vstddata)), ncol = 1) - c
+ if (method == "OLS") {
+ Z = try(solve(t(X) %*% X), silent = T)
+ if (inherits(Z, "try-error")) {
+ print("X'X is not invertible. Switch to TML")
+ }
+ else {
+ theta = solve(t(X) %*% X) %*% t(X) %*% vy
+ rm(X)
+ rm(vy)
+ }
+ }
+ if (method == "TML") {
+ inittheta = rep(0, dim(X)[2])
+ l = -2.272
+ u = 1.6675
+ nonoutliers = c(1:length(vy))[(firststepresids >
+ l) & (firststepresids < u)]
+ truncvy = vy[nonoutliers]
+ rm(vy)
+ truncX = X[nonoutliers, ]
+ rm(X)
+ negtruncLLH = function(theta) {
+ res = truncvy - truncX %*% matrix(theta, ncol = 1)
+ return(mean(-res - c + exp(2 * (res + c))/2))
+ }
+ grnegtruncLLH = function(theta) {
+ res = truncvy - truncX %*% matrix(theta, ncol = 1)
+ dres = -truncX
+ return(apply(-dres + as.vector(exp(2 * (res +
+ c))) * dres, 2, "mean"))
+ }
+ est = optim(par = inittheta, fn = negtruncLLH, gr = grnegtruncLLH,
+ method = "BFGS")
+ theta = est$par
+ rm(truncX)
+ rm(truncvy)
+ }
+ plot(seas, main = "Non-parametric (dashed line) and parametric (full line) periodicity",
+ xlab = "intraday period", type = "l", lty = 3)
+ seas = RTAQ:::diurnalfit(theta = theta, P1 = P1, P2 = P2, intraT = intraT,
+ dummies = dummies)
+ lines(seas, lty = 1)
+ return(list(seas, theta))
+ }
+ else {
+ return(list(seas))
+ }
}
diurnalfit = function( theta , P1 , P2 , intraT , dummies=F )
Modified: pkg/RTAQ/man/spotVol.rd
===================================================================
--- pkg/RTAQ/man/spotVol.rd 2011-03-12 18:52:40 UTC (rev 573)
+++ pkg/RTAQ/man/spotVol.rd 2011-03-13 18:45:47 UTC (rev 574)
@@ -22,20 +22,12 @@
robust regression. The regression specification consists either of one dummy for each intraday period (dummies=TRUE) or the flexible fourrier
form with P1 cosinus and P2 sinus terms. For more details on the classical methods, see Taylor and Xu (1997) and Andersen et al. (1997).
For the jump robust versions, see Boudt et al. (2010).
-
}
\usage{
-spotVol(pdata,
- dailyvol = "bipower",
- periodicvol = "TML",
- on = "minutes",
- k = 5,
- dummies = FALSE ,
- P1 = 4,
- P2 = 2,
- ...)
-}
+spotVol(pdata, dailyvol = "bipower", periodicvol = "TML",
+ on = "minutes", k = 5, dummies = FALSE, P1 = 4, P2 = 2,
+ marketopen = "09:30:00", marketclose = "16:00:00")}
\arguments{
\item{pdata}{xts object, containing the price series.}
@@ -50,7 +42,8 @@
If it false, the parametric estimator uses the Flexible Fourrier specification. FALSE by default.}
\item{P1}{ is a positive integer valued parameter that corresponds to the number of cosinus terms used in the flexible fourrier specification for the periodicity function, see Andersen et al. (1997) for details.}
\item{P2}{ is a positive integer valued parameter that corresponds to the number of sinus terms used in the flexible fourrier specification for the periodicity function, see Andersen et al. (1997) for details.}
- \item{...}{additional arguments}
+ \item{marketopen}{the market opening time, by default: marketopen = "09:30:00".}
+ \item{marketclose}{the market closing time, by default: marketclose = "16:00:00".}
}
\details{
@@ -82,11 +75,10 @@
data("sample_real5minprices");
#compute and plot intraday periodicity
-out = spotVol(price,P1=6,P2=4,method="OLS",k=5, dummies=FALSE);
+out = spotVol(price,P1=6,P2=4,periodicvol="TML",k=5, dummies=FALSE);
head(out);
}
-
\keyword{ volatility}
\author{ Jonathan Cornelissen and Kris Boudt}
More information about the Blotter-commits
mailing list