[Blotter-commits] r1032 - in pkg/RTAQ: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue May 29 02:58:55 CEST 2012
Author: jonathan
Date: 2012-05-29 02:58:54 +0200 (Tue, 29 May 2012)
New Revision: 1032
Removed:
pkg/RTAQ/R/realized.R
pkg/RTAQ/man/as.realizedObject.Rd
Modified:
pkg/RTAQ/DESCRIPTION
pkg/RTAQ/NAMESPACE
pkg/RTAQ/R/HAR_model.R
pkg/RTAQ/man/harModel.Rd
Log:
updated harModel, jump aggregation and leverage effects can now be included!
Modified: pkg/RTAQ/DESCRIPTION
===================================================================
--- pkg/RTAQ/DESCRIPTION 2012-05-28 21:32:03 UTC (rev 1031)
+++ pkg/RTAQ/DESCRIPTION 2012-05-29 00:58:54 UTC (rev 1032)
@@ -1,7 +1,7 @@
Package: RTAQ
Type: Package
Title: RTAQ: Tools for the analysis of trades and quotes in R
-Version: 0.2
+Version: 0.3
Author: Jonathan Cornelissen, Kris Boudt
Maintainer: Jonathan Cornelissen <Jonathan.cornelissen at econ.kuleuven.be>
Description: The Trades and Quotes data of the New York Stock Exchange is a popular input for the implementation of intraday trading strategies, the measurement of liquidity and volatility and investigation of the market microstructure, among others. This package contains a collection of R functions to carefully clean and match the trades and quotes data, calculate ex post liquidity and volatility measures and detect price jumps in the data.
Modified: pkg/RTAQ/NAMESPACE
===================================================================
--- pkg/RTAQ/NAMESPACE 2012-05-28 21:32:03 UTC (rev 1031)
+++ pkg/RTAQ/NAMESPACE 2012-05-29 00:58:54 UTC (rev 1032)
@@ -28,7 +28,6 @@
aggregateQuotes,
aggregateTrades,
aggregatets,
-as.realizedObject,
matchTradesQuotes,
TAQLoad,
refreshTime,
Modified: pkg/RTAQ/R/HAR_model.R
===================================================================
--- pkg/RTAQ/R/HAR_model.R 2012-05-28 21:32:03 UTC (rev 1031)
+++ pkg/RTAQ/R/HAR_model.R 2012-05-29 00:58:54 UTC (rev 1032)
@@ -30,7 +30,8 @@
return(zstat);
}
- harModel = function(data, periods = c(1,5,22), RVest = c("RCov","RBPCov"), type="HARRV", jumptest="ABDJumptest",alpha=0.05,h=1,transform=NULL,...){
+ harModel = function(data, periods = c(1,5,22), periodsJ = c(1,5,22), leverage=NULL, RVest = c("RCov","RBPCov"), type="HARRV",
+ jumptest="ABDJumptest",alpha=0.05,h=1,transform=NULL, ...){
nperiods = length(periods); # Number of periods to aggregate over
nest = length(RVest); # Number of RV estimators
if( !is.null(transform) ){ Ftransform = match.fun(transform); }
@@ -40,95 +41,100 @@
# Get the daily RMs (in a non-robust and robust way)
RV1 = match.fun( RVest[1]);
RM1 = apply.daily( data, RV1 );
- #save dates:
+ # save dates:
alldates = index(RM1)
if( nest == 2 ){
RV2 = match.fun( RVest[2]);
RM2 = apply.daily( data, RV2 ); }
}
-
-
+
if( sum(data<0) == 0 ){ #The input is most likely already realized measures
- dimdata = dim(data)[2];
- alldates = index(data);
- RM1 = data[,1];
- if( dimdata > 1 ){ RM2 = data[,2]; }
- if( type != "HARRV" ){ warning("Please provide returns as input for the type of model you want to estimate. All your returns are positive which is quite unlikely honestly. Only for the HAR-RV model you can input realized measures.") }
+ dimdata = dim(data)[2];
+ alldates = index(data);
+ RM1 = data[,1];
+ if( dimdata > 1 ){ RM2 = data[,2]; }
+ if( type != "HARRV" ){ warning("Please provide returns as input for the type of model you want to estimate. All your returns are positive which is quite unlikely honestly. Only for the HAR-RV model you can input realized measures.") }
}
# Get the matrix for estimation of linear model
- maxp = max(periods); #Number of aggregation levels
+ maxp = max(periods,periodsJ); #max number of aggregation levels
+ if(!is.null(leverage)){ maxp = max(maxp,leverage) }
n = length(RM1); #Number of Days
- RVmatrix1 = matrix(nrow=n,ncol=nperiods);
-
- for(i in 1:nperiods){
- if(periods[i]==1){ RVmatrix1[,i] = RM1
- }else{ RVmatrix1[(periods[i]:n),i] = rollmean(x=RM1,k=periods[i],align="left") }
- } #end loop over periods for standard RV estimator
- colnames(RVmatrix1) = paste("RV",periods,sep="");
- # Aggregate and subselect y
- if( h == 1 ){ y = RM1[(maxp+1):n]; }
- if( h != 1 ){
- y = matrix( nrow=length(RM1), ncol=1 ); colnames(y) = "y";
- y[(h:n),] = rollmean(x=RM1,k=h,align="left");
- y = matrix(y[((maxp+h):n),],ncol=1); y=as.data.frame(y) }
+ # Aggregate RV:
+ RVmatrix1 = aggRV(RM1,periods);
+ if( nest==2 ){ RVmatrix2 = aggRV(RM2,periods); } # In case a jumprobust estimator is supplied
- # Only keep useful parts:
- x1 = RVmatrix1[(maxp:(n-h)),];
+ # Aggregate and subselect y:
+ y = aggY(RM1,h,maxp);
- # TODO: add transformations here (srqr,log,..) see paper
- if(nest==2){ # In case a jumprobust estimator is supplied
- RVmatrix2 = matrix(nrow=n,ncol=nperiods);
- for(i in 1:nperiods){
- if(periods[i]==1){ RVmatrix2[,i] = RM2;
- }else{ RVmatrix2[(periods[i]:n),i] = rollmean(x=RM2,k=periods[i],align="left") }
- colnames(RVmatrix2) = paste("RV",periods,sep="");
- x2 = RVmatrix2[(maxp:(n-h)),];
- } #end loop over periods for robust RV estimator
- }
-
- # Estimate the model parameters, according to type of model :
- # First model type: traditional HAR-RV:
- if( type == "HARRV" ){
- if(!is.null(transform)){ y = Ftransform(y); x1 = Ftransform(x1) }
+ # Only keep useful parts:
+ x1 = RVmatrix1[(maxp:(n-h)),];
+ if( nest==2 ){ x2 = RVmatrix2[(maxp:(n-h)),]; } # In case a jumprobust estimator is supplied
+
+ # Jumps:
+ if(type!="HARRV"){ # If model type is as such that you need jump component
+ J = pmax( RM1 - RM2,0 ); # Jump contributions should be positive
+ J = aggJ(J,periodsJ);
+ }
+
+ if( !is.null(leverage) ){
+ if( sum(data<0) == 0 ){ warning("You cannot use leverage variables in the model in case your input consists of Realized Measures") }
+ # Get close-to-close returns
+ e = apply.daily(data,sum); #Sum logreturns daily
+ # Get the rmins:
+ rmintemp = pmin(e,0);
+ # Aggregate everything:
+ rmin = aggRV(rmintemp,periods=leverage,type="Rmin");
+ # Select:
+ rmin = rmin[(maxp:(n-h)),];
+ }else{ rmin = matrix(ncol=0,nrow=dim(x1)[1]) }
+
+ ###############################
+ # Estimate the model parameters, according to type of model :
+ # First model type: traditional HAR-RV:
+ if( type == "HARRV" ){
+ if(!is.null(transform)){ y = Ftransform(y); x1 = Ftransform(x1) }
+ x1 = cbind(x1,rmin);
model = estimhar(y=y,x=x1);
model$transform = transform; model$h = h; model$type = "HARRV"; model$dates = alldates[(maxp+h):n];
class(model) = c("harModel","lm");
return( model )
- } #End HAR-RV if cond
+ } #End HAR-RV if cond
if( type == "HARRVJ" ){
- J = pmax( RM1 - RM2,0 ); #Jump contributions should be positive
- J = matrix(J[(maxp:(n-h)),]); colnames(J) = "J";
+ J = J[(maxp:(n-h)),];
x = cbind(x1,J); # bind jumps to RV data
if(!is.null(transform)){ y = Ftransform(y); x = Ftransform(x); }
- model = estimhar(y=y,x=x);
+ x = cbind(x,rmin);
+ model = estimhar(y=y,x=x);
model$transform = transform; model$h = h; model$type = "HARRVJ"; model$dates = alldates[(maxp+h):n];
class(model) = c("harModel","lm");
return( model )
}#End HAR-RV-J if cond
- if( type == "HARRVCJ" ){
- # Get the jumps:
- J = pmax( RM1 - RM2,0 ); # Jump contributions should be positive
+ if( type == "HARRVCJ" ){
# Are the jumps significant? if not set to zero:
if( jumptest=="ABDJumptest" ){
-
TQ = apply.daily(data, TQfun);
+ J = J[,1];
teststats = ABDJumptest(RV=RM1,BPV=RM2,TQ=TQ );
}else{ jtest = match.fun(jumptest); teststats = jtest(data,...) }
Jindicators = teststats > qnorm(1-alpha);
J[!Jindicators] = 0;
- J = matrix(J[(maxp:(n-h)),]); colnames(J) = "J"; Jindicators = Jindicators[(maxp:(n-h))];
- # Get continuus components if necessary RV measures if necessary:
- Cmatrix = matrix( nrow = dim(x1)[1], ncol = dim(x1)[2] );
- Cmatrix[Jindicators,] = x2[Jindicators,]; #Fill with robust one in case of jump
- Cmatrix[(!Jindicators),] = x1[(!Jindicators),]; #Fill with non-robust one in case of no-jump
- colnames(Cmatrix) = paste("C",periods,sep="");
-
- x = cbind(Cmatrix,J); # bind jumps to RV data
- if(!is.null(transform)){ y = Ftransform(y); x = Ftransform(x); }
+ # Get continuus components if necessary RV measures if necessary:
+ Cmatrix = matrix( nrow = dim(RVmatrix1)[1], ncol = 1 );
+ Cmatrix[Jindicators,] = RVmatrix2[Jindicators,1]; #Fill with robust one in case of jump
+ Cmatrix[(!Jindicators)] = RVmatrix1[(!Jindicators),1]; #Fill with non-robust one in case of no-jump
+ # Aggregate again:
+ Cmatrix <- aggRV(Cmatrix,periods,type="C");
+ Jmatrix <- aggJ(J,periodsJ);
+ # subset again:
+ Cmatrix <- Cmatrix[(maxp:(n-h)),];
+ Jmatrix <- Jmatrix[(maxp:(n-h)),];
+ x = cbind(Cmatrix,Jmatrix); # bind jumps to RV data
+ if(!is.null(transform)){ y = Ftransform(y); x = Ftransform(x); }
+ x = cbind(x,rmin);
model = estimhar( y=y, x=x );
model$transform = transform; model$h = h; model$type = "HARRVCJ"; model$dates = alldates[(maxp+h):n];
class(model) = c("harModel","lm");
@@ -136,7 +142,7 @@
}
} #End function harModel
-
+ #################################################################
estimhar = function(y, x){ #Potentially add stuff here
colnames(y)="y";
output = lm( formula(y~x), data=cbind(y,x));
@@ -158,6 +164,42 @@
return(list(modeldescription,betas))
}
+ aggRV <- function(RM1,periods,type="RV"){
+ n = length(RM1);
+ nperiods = length(periods);
+ RVmatrix1 = matrix(nrow=n,ncol=nperiods);
+ for(i in 1:nperiods){
+ if(periods[i]==1){ RVmatrix1[,i] = RM1;
+ }else{ RVmatrix1[(periods[i]:n),i] = rollmean(x=RM1,k=periods[i],align="left") }
+ } #end loop over periods for standard RV estimator
+ colnames(RVmatrix1) = paste(type,periods,sep="");
+ return(RVmatrix1);
+ }
+
+ aggJ <- function( J, periodsJ ){
+ n = length(J);
+ nperiods = length(periodsJ);
+ JM = matrix(nrow=n,ncol=nperiods);
+ for(i in 1:nperiods){
+ if(periodsJ[i]==1){ JM[,i] = J;
+ }else{ JM[(periodsJ[i]:n),i] = rollmean( x=J, k=periodsJ[i], align="left") }
+ } # End loop over periods for standard RV estimator
+ colnames(JM) = paste("J",periodsJ,sep="");
+ return(JM)
+ }
+
+ aggY = function(RM1,h,maxp){
+ n = length(RM1);
+ if( h == 1 ){ y = RM1[(maxp+1):n]; }
+ if( h != 1 ){
+ y = matrix( nrow=length(RM1), ncol=1 ); colnames(y) = "y";
+ y[(h:n),] = rollmean(x=RM1,k=h,align="left");
+ y = matrix(y[((maxp+h):n),],ncol=1); y=as.data.frame(y) }
+ return(y);
+ }
+
+
+#########################################################################
# Print method for harmodel:
print.harModel = function(x, digits = max(3, getOption("digits") - 3), ...){
formula = getHarmodelformula(x); modeldescription = formula[[1]]; betas = formula[[2]];
@@ -202,6 +244,7 @@
observed = x$model$y;
fitted = x$fitted.values;
dates = x$dates;
+ dates = as.POSIXct(dates);
observed = xts(observed, order.by=dates);
fitted = xts(fitted, order.by=dates);
type = x$type;
@@ -214,7 +257,6 @@
# axis(1,time(b)[ind], format(time(b)[ind],), las=2, cex.axis=0.8); not used anymore
# axis(2);
lines(fitted,col="blue",lwd=2);
- legend("topleft", c("Observed RV","Forecasted RV"), cex=0.8, col=c("red","blue"),lty=1, lwd=2, bty="n");
+ legend("topleft", c("Observed RV","Forecasted RV"), cex=1.1, col=c("red","blue"),lty=1, lwd=2, bty="n");
}
-
\ No newline at end of file
Deleted: pkg/RTAQ/R/realized.R
===================================================================
--- pkg/RTAQ/R/realized.R 2012-05-28 21:32:03 UTC (rev 1031)
+++ pkg/RTAQ/R/realized.R 2012-05-29 00:58:54 UTC (rev 1032)
@@ -1,45 +0,0 @@
-to.milliseconds <- function(td, startDate) {
- # function to extract milliseconds from “timeDate” time index
- # Based on initial code from Eric Zivot
-
- ms = as.numeric(td at Data) - as.numeric(startDate at Data)
- return(ms*1000)
-}
-
-
-as.realizedObject = function( ts, cts = FALSE , makeReturns = TRUE ){
-require('realized');
- # takes the RTAQ data as input and outputs a list with each element a realizedobj of the corresponding day
- # Based on initial code from Eric Zivot
- days = unique( format( time(ts) , "%Y-%m-%d"))
- i = 1;
- out = list();
-
- #one day only
- if(length(days)==1){
- startDate = as.timeDate(paste( days, "00:00:00")); #"%Y-%m-%d %H:%M:%S"
- temp = list( data = as.character(ts[days]) ,
- milliseconds = to.milliseconds(td=index(ts[days]), startDate) );
- out = realizedObject( x = temp , makeReturns=makeReturns, cts=cts,
- millisstart=9.5*60*60*1000, millisend=16*60*60*1000)
- }
-
- #multiple days: make list
- if(length(days)!=1){
- for( day in days ){
-
- startDate = as.timeDate(paste( day, "00:00:00")); #"%Y-%m-%d %H:%M:%S"
-
- temp = list( data = as.character(ts[day]) ,
- milliseconds = to.milliseconds(td=index(ts[day]), startDate) )
-
- out[[i]] = realizedObject( x = temp , makeReturns=makeReturns, cts=cts,
- millisstart=9.5*60*60*1000, millisend=16*60*60*1000)
- # important to set cts to false, otherwise all zero prices!
- i = i + 1;
- }
- }
-
- #if( cts ){ print("cts=TRUE: for the milliseconds without observations prices are inserted")}
- return(out)
-}
\ No newline at end of file
Deleted: pkg/RTAQ/man/as.realizedObject.Rd
===================================================================
--- pkg/RTAQ/man/as.realizedObject.Rd 2012-05-28 21:32:03 UTC (rev 1031)
+++ pkg/RTAQ/man/as.realizedObject.Rd 2012-05-29 00:58:54 UTC (rev 1032)
@@ -1,34 +0,0 @@
-\name{as.realizedObject}
-\Rdversion{1.1}
-\alias{as.realizedObject}
-\title{Convert xts object to a realized object}
-
-\description{Function converts an xts timeseries (one dimensional, possibly multiple days) into
-a realized object from the package "realized" by Scott Payseur.}
-
-\usage{as.realizedObject( ts, cts = FALSE , makeReturns = TRUE )}
-
-\arguments{
-\item{ts}{ xts object, containing the original time series (one dimensional, possibly multiple days).}
-\item{cts}{ Calendar time sampling or tick time sampling. See realizedObject function in "realized" package by Scott Payseur.}
-\item{makeReturns}{ boolean, indicates whether you pass prices. TRUE by default. See realizedObject function in "realized" package by Scott Payseur.}
-}
-
-\value{A realized object.}
-
-\author{ Jonathan Cornelissen and Kris Boudt. Initial code by Eric Zivot}
-\keyword{data manipulation}
-
-\examples{
-#one day of data:
-data("sample_tdata");
-x = as.realizedObject(sample_tdata$PRICE);
-
-#multiple days:
-data("sample_5minprices");
-ts1 = sample_5minprices$stock1;
-#construct realized object:
-x = as.realizedObject(ts1);
-#x now contains a list with a realized object for each day
-#see manual of the "realized" package by Scott Payseur for more info.
-}
\ No newline at end of file
Modified: pkg/RTAQ/man/harModel.Rd
===================================================================
--- pkg/RTAQ/man/harModel.Rd 2012-05-28 21:32:03 UTC (rev 1031)
+++ pkg/RTAQ/man/harModel.Rd 2012-05-29 00:58:54 UTC (rev 1032)
@@ -10,14 +10,18 @@
\usage{
- harModel(data, periods = c(1, 5, 22), RVest = c("RCov", "RBPCov"),
+ harModel(data, periods = c(1, 5, 22), periodsJ = c(1,5,22), leverage=NULL, RVest = c("RCov", "RBPCov"),
type = "HARRV", jumptest = "ABDJumptest", alpha = 0.05, h = 1,
transform = NULL, ...) }
\arguments{
- \item{data}{ an xts-object containing the intraday returns.}
- \item{periods}{ a vector of integers indicating over how days the realized measures in the model should be aggregated. By default periods = c(1,5,22), which corresponds to one day, one week and one month respectively. This default is in line with Andersen et al. (2007).}
- \item{ RVest}{ a character vector with one or two elements.
+ \item{data}{ an xts-object containing the intraday (log-)returns.}
+ \item{periods}{ a vector of integers indicating over how days the realized measures in the model should be aggregated. By default periods = c(1,5,22), which corresponds to one day, one week and one month respectively. This default is in line with Andersen et al. (2007).}
+ \item{periodsJ}{ a vector of integers indicating over what time periods the jump components in the model should be aggregated. By default periodsJ = c(1,5,22), which corresponds to one day, one week and one month respectively.}
+ \item{leverage}{ a vector of integers indicating over what periods the negative returns should be aggregated.
+ See Corsi and Reno (2012) for more information. By default leverage=NULL and the model assumes the absence of a leverage effect. Set leverage= c(1,5,22) to mimic the analysis in Corsi and Reno (2012).
+ }
+ \item{RVest}{ a character vector with one or two elements.
The first element refers to the name of the function to estimate the daily integrated variance (non-jump-robust), while the second element refers to the name of the function to estimate the continuous component of daily volatility (jump-robust). By default RVest = c("RCov","RBPCov"), i.e. using the Realized Volatility and Realized Bi-Power Variance.}
\item{type}{ a string referring to the type of HAR model you would like to estimate. By default type = "HARRV", the most basic model. Other valid options are type = "HARRVJ" or type = "HARRVCJ".}
\item{jumptest}{ the function name of a function used to test whether the test statistic which determines whether the jump variability is significant that day. By default jumptest = "ABDJumptest", hence using the test statistic in Equation or Equation (18) of Andersen et al. (2007).}
@@ -36,6 +40,19 @@
The function outputs an object of class \code{harModel} and \code{\link{lm}} (so \code{harModel} is a subclass of \code{\link{lm}}). So far I only added a print method as you can see in the examples. Input here is welcome, what should a plot of an "harmodel" object look like? What other methods are useful?
}
+
+\references{
+Andersen, T. G., T. Bollerslev, and F. Diebold (2007). Roughing it up: includ-
+ing jump components in the measurement, modelling and forecasting of return
+volatility. The Review of Economics and Statistics 89, 701-720.
+
+Corsi, F. (2009). A simple approximate long memory model of realized volatility.
+Journal of Financial Econometrics 7, 174-196.
+
+Corsi, F. and Reno R. (2012). Discrete-time volatility forecasting with persistent leverage effect and the link with continuous-time volatility modeling. Journal of Business and Economic Statistics, forthcoming.
+}
+
+
\author{ Jonathan Cornelissen and Kris Boudt}
\keyword{forecasting}
@@ -44,11 +61,11 @@
data("sample_5minprices_jumps");
data = sample_5minprices_jumps[,1];
data = makeReturns(data); #Get the high-frequency return data
-
- x = harModel(data, periods = c(1,5,10), RVest = c("RCov","RBPCov"),
- type="HARRVCJ",h=3,transform="sqrt"); # Estimate the HAR model of type HARRVCJ
+
+ x = harModel(data, periods = c(1,5,10), periodsJ=c(1,5,10), RVest = c("RCov","RBPCov"),
+ type="HARRVCJ",transform="sqrt"); # Estimate the HAR model of type HARRVCJ
class(x);
- x
+ x
##### Example 2: #####
# Forecasting daily Realized volatility for DJI 2008 using the basic harModel: HARRV
More information about the Blotter-commits
mailing list