[Blotter-commits] r974 - pkg/RTAQ/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Mar 15 19:03:42 CET 2012
Author: kboudt
Date: 2012-03-15 19:03:42 +0100 (Thu, 15 Mar 2012)
New Revision: 974
Modified:
pkg/RTAQ/R/volatility.R
Log:
Update two time scale functions, whereby different time scales can be used for different assets. Also corrected bug reported by Darko.
Modified: pkg/RTAQ/R/volatility.R
===================================================================
--- pkg/RTAQ/R/volatility.R 2012-03-15 17:11:14 UTC (rev 973)
+++ pkg/RTAQ/R/volatility.R 2012-03-15 18:03:42 UTC (rev 974)
@@ -300,7 +300,8 @@
# return(rcor);
#}
-RTSRV = function (pdata, startIV = NULL, noisevar = NULL, K = 300, J = 1, eta = 9){
+RTSRV = function (pdata, startIV = NULL, noisevar = NULL, K = 300, J = 1,
+ eta = 9){
logprices = log(as.numeric(pdata))
n = length(logprices)
nbarK = (n - K + 1)/(K)
@@ -315,31 +316,33 @@
logreturns_K = c(logreturns_K, diff(logprices[sel]))
vdelta_K = c(vdelta_K, diff(seconds[sel])/secday)
}
-
for (j in 1:J) {
sel = seq(j, n, J)
logreturns_J = c(logreturns_J, diff(logprices[sel]))
vdelta_J = c(vdelta_J, diff(seconds[sel])/secday)
}
if (is.null(noisevar)) {
- logreturns_1 = diff(logprices)
- noisevar = 1/(2 * n) * (sum(logreturns_1^2) - RTAQ:::TSRV(pdata,
- K = K, J = J))
+ noisevar = max(0,1/(2 * nbarJ) * (sum(logreturns_J^2)/J - RTAQ:::TSRV(pdata1,K=K,J=J)))
}
-
if (!is.null(startIV)) {
RTSRV = startIV
}
if (is.null(startIV)) {
- RTSRV = (1/K) * MedRV(logreturns_K)
+ sel = seq(1, n, K)
+ RTSRV = MedRV(diff(logprices[sel]))
}
-
iter = 1
while (iter <= 20) {
- I_K = 1 * (logreturns_K^2 <= eta * (RTSRV * vdelta_K + 2 * noisevar))
- I_J = 1 * (logreturns_J^2 <= eta * (RTSRV * vdelta_J + 2 * noisevar))
- if( sum(I_J)==0 ){I_J = rep(1,length(logreturns_J));}
- if( sum(I_K)==0 ){I_K = rep(1,length(logreturns_K));}
+ I_K = 1 * (logreturns_K^2 <= eta * (RTSRV * vdelta_K +
+ 2 * noisevar))
+ I_J = 1 * (logreturns_J^2 <= eta * (RTSRV * vdelta_J +
+ 2 * noisevar))
+ if (sum(I_J) == 0) {
+ I_J = rep(1, length(logreturns_J))
+ }
+ if (sum(I_K) == 0) {
+ I_K = rep(1, length(logreturns_K))
+ }
RTSRV = adj * (zeta * (1/K) * sum(logreturns_K^2 * I_K)/mean(I_K) -
((nbarK/nbarJ) * zeta * (1/J) * sum(logreturns_J^2 *
I_J)/mean(I_J)))
@@ -372,157 +375,231 @@
-RTSCov = function (pdata,cor=FALSE, startIV=NULL, noisevar = NULL, K = 300, J = 1, eta = 9, makePsd=FALSE)
-{
- #pdata is list with each item an xts object with 1 tick-by-tick timeseries
- # if(hasArg(data)){ rdata = data }
- # if(makeReturns){ rdata = makeReturns(rdata)}
- if(!is.list(pdata)){ n = 1 }else{
- n = length(pdata);
- if(n==1){pdata=pdata[[1]];}
- }
-
- if( n == 1 ){ return( RTSRV( pdata, startIV=startIV , noisevar = noisevar , K = K , J = J, eta = eta ))}
-
- if( n > 1 ){
-
- cov = matrix(rep(0, n * n), ncol = n)
- diagonal = c()
- for (i in 1:n){
- diagonal[i] = RTSRV(pdata[[i]], startIV = startIV[i], noisevar = noisevar[i] , K = K , J = J , eta = eta)
- }
- diag(cov) = diagonal
- for (i in 2:n) {
- for (j in 1:(i - 1)) {
- cov[i, j] = cov[j, i] = RTSCov_bi(pdata[[i]], pdata[[j]],
- startIV1 = diagonal[i], startIV2 = diagonal[j],
- noisevar1 = noisevar[i], noisevar2 = noisevar[j],
- K = K, J = J, eta = eta)
- }
- }
-
- if(cor==FALSE){
- if(makePsd==TRUE){cov = makePsd(cov);}
- return(cov)}
- if(cor==TRUE){
- invsdmatrix = try(solve(sqrt(diag(diag(cov)))),silent=F);
- if( inherits( invsdmatrix, "try-error") ){
- rcor = invsdmatrix%*%cov%*%invsdmatrix;
- if(makePsd==TRUE){rcor = makePsd(rcor);}
- return(rcor)}
- }
- }
+RTSCov = function (pdata, cor = FALSE, startIV = NULL, noisevar = NULL,
+ K = 300, J = 1,
+ K_cov = NULL , J_cov = NULL,
+ K_var = NULL , J_var = NULL ,
+ eta = 9, makePsd = FALSE){
+ if (!is.list(pdata)) {
+ n = 1
+ }
+ else {
+ n = length(pdata)
+ if (n == 1) {
+ pdata = pdata[[1]]
+ }
+ }
+ if (n == 1) {
+ return(RTAQ:::RTSRV(pdata, startIV = startIV, noisevar = noisevar,
+ K = K, J = J, eta = eta))
+ }
+ if (n > 1) {
+ cov = matrix(rep(0, n * n), ncol = n)
+ diagonal = c()
+ if( is.null(K_cov)){ K_cov = K }
+ if( is.null(J_cov)){ J_cov = J }
+ if( is.null(K_var)){ K_var = rep(K,n) }
+ if( is.null(J_var)){ J_var = rep(J_var,n) }
+ for (i in 1:n) {
+ diagonal[i] = RTAQ:::RTSRV(pdata[[i]], startIV = startIV[i],
+ noisevar = noisevar[i], K = K_var[i], J = J_var[i],
+ K_cov = K_cov , J_cov = J_cov , eta = eta)
+ }
+ diag(cov) = diagonal
+ if( is.null(K_cov)){ K_cov = K }
+ if( is.null(J_cov)){ J_cov = J }
+ for (i in 2:n) {
+ for (j in 1:(i - 1)) {
+ cov[i, j] = cov[j, i] = RTSCov_bi(pdata[[i]],
+ pdata[[j]], startIV1 = diagonal[i], startIV2 = diagonal[j],
+ noisevar1 = noisevar[i], noisevar2 = noisevar[j],
+ K = K, J = J, eta = eta)
+ }
+ }
+ if (cor == FALSE) {
+ if (makePsd == TRUE) {
+ cov = makePsd(cov)
+ }
+ return(cov)
+ }
+ if (cor == TRUE) {
+ invsdmatrix = try(solve(sqrt(diag(diag(cov)))), silent = F)
+ if (!inherits(invsdmatrix, "try-error")) {
+ rcor = invsdmatrix %*% cov %*% invsdmatrix
+ if (makePsd == TRUE) {
+ rcor = makePsd(rcor)
+ }
+ return(rcor)
+ }
+ }
+ }
}
-RTSCov_bi = function(pdata1,pdata2,startIV1=NULL,startIV2=NULL, noisevar1 = NULL, noisevar2 = NULL, K = 300, J = 1, eta = 9){
- # pdata is an xts object containing tick-by-tick prices:
-
- #get refresh prices:
- x = refreshTime(list(pdata1,pdata2));
- newprice1 = x[,1];
- newprice2 = x[,2];
+RTSCov_bi =
+function (pdata1, pdata2, startIV1 = NULL, startIV2 = NULL, noisevar1 = NULL,
+ noisevar2 = NULL, K = 300, J = 1,
+ K_cov = NULL , J_cov = NULL ,
+ K_var1 = NULL , K_var2 = NULL ,
+ J_var1 = NULL , J_var2 = NULL ,
+ eta = 9)
+{
+ if( is.null(K_cov)){ K_cov = K } ; if( is.null(J_cov)){ J_cov = J }
+ if( is.null(K_var1)){ K_var1 = K } ; if( is.null(K_var2)){ K_var2 = K }
+ if( is.null(J_var1)){ J_var1 = J } ; if( is.null(J_var2)){ J_var2 = J }
- #get the rest:
- logprices1 = log(as.numeric(newprice1));
- logprices2 = log(as.numeric(newprice2));
- seconds = as.numeric(as.POSIXct(index(newprice1)));
+ # Calculation of the noise variance and TSRV for the truncation
+
- n = length(logprices1);
- nbarK = (n - K + 1)/(K);
- nbarJ = (n - J + 1)/(J);
- adj = n/((K-J) * nbarK);
- logreturns_K1 = logreturns_K2 = vdelta_K = logreturns_J1 = logreturns_J2 = vdelta_J = c();
- secday = last(seconds) - first(seconds);
- for( k in 1:K ){
- sel.avg = seq(k,n,K)
- logreturns_K1 = c( logreturns_K1 , diff( logprices1[sel.avg]));
- logreturns_K2 = c( logreturns_K2 , diff( logprices2[sel.avg]));
- vdelta_K = c( vdelta_K , diff( seconds[sel.avg])/secday)
- }
+ if ( is.null(noisevar1) ) {
+ logprices1 = log(as.numeric(pdata1))
+ n_var1 = length(logprices1)
+ nbarK_var1 = (n_var1 - K_var1 + 1)/(K_var1) ;
+ nbarJ_var1 = (n_var1 - J_var1 + 1)/(J_var1)
+ adj_var1 = n_var1/((K_var1 - J_var1) * nbarK_var1)
+
+ logreturns_K1 = logreturns_J1 = c()
+ for (k in 1:K_var1) {
+ sel.avg = seq(k, n_var1, K_var1)
+ logreturns_K1 = c(logreturns_K1, diff(logprices1[sel.avg]))
+ }
+ for (j in 1:J_var1) {
+ sel.avg = seq(j, n_var1, J_var1)
+ logreturns_J1 = c(logreturns_J1, diff(logprices1[sel.avg]))
+ }
+ if( is.null(noisevar1) ){
+ noisevar1 = max(0,1/(2 * nbarJ_var1) * (sum(logreturns_J1^2)/J_var1 - RTAQ:::TSRV(pdata1,K=K_var1,J=J_var1)))
+ }
+ }
+ if (is.null(noisevar2)) {
+ logprices2 = log(as.numeric(pdata2))
+ n_var2 = length(logprices2)
+ nbarK_var2 = (n_var2 - K_var2 + 1)/(K_var2) ;
+ nbarJ_var2 = (n_var2 - J_var2 + 1)/(J_var2)
+ adj_var2 = n_var2/((K_var2 - J_var2) * nbarK_var2)
+
+ logreturns_K2 = logreturns_J2 = c()
+ for (k in 1:K_var2) {
+ sel.avg = seq(k, n_var2, K_var2)
+ logreturns_K2 = c(logreturns_K2, diff(logprices2[sel.avg]))
+ }
+ for (j in 1:J_var2) {
+ sel.avg = seq(j, n_var2, J_var2)
+ logreturns_J2 = c(logreturns_J2, diff(logprices2[sel.avg]))
+ }
+ noisevar2 = max(0,1/(2 * nbarJ_var2) * (sum(logreturns_J2^2)/J_var2 - RTAQ:::TSRV(pdata2,K=K_var2,J=J_var2)))
+ }
- for( j in 1:J ){
- sel.avg = seq(j,n,J)
- logreturns_J1 = c( logreturns_J1 , diff( logprices1[sel.avg]));
- logreturns_J2 = c( logreturns_J2 , diff( logprices2[sel.avg]));
- vdelta_J = c( vdelta_J , diff( seconds[sel.avg])/secday)
- }
-
- if( is.null(noisevar1) ){
- logreturns1 = diff( logprices1 );
- noisevar1 = 1/(2*n)*( sum( logreturns1^2 ) - TSRV(pdata1) );
- }
+ if (!is.null(startIV1)) {
+ RTSRV1 = startIV1
+ }else{
+ RTSRV1 = RTSRV(pdata1, noisevar = noisevar1, K = K_var1, J = J_var1, eta = eta)
+ }
+ if (!is.null(startIV2)) {
+ RTSRV2 = startIV2
+ }else{
+ RTSRV2 = RTSRV(pdata2, noisevar = noisevar2, K = K_var2, J = J_var2, eta = eta)
+ }
+
+ # Refresh time is for the covariance calculation
+
+ x = refreshTime(list(pdata1, pdata2))
+ newprice1 = x[, 1]
+ newprice2 = x[, 2]
+ logprices1 = log(as.numeric(newprice1))
+ logprices2 = log(as.numeric(newprice2))
+ seconds = as.numeric(as.POSIXct(index(newprice1)))
+ secday = last(seconds) - first(seconds)
+ K = K_cov ; J = J_cov ;
- if( is.null(noisevar2) ){
- logreturns2 = diff( logprices2 );
- noisevar2 = 1/(2*n)*( sum( logreturns2^2 ) - TSRV(pdata2) );
- }
-
- zeta = 1/pchisq( eta , 3 )
-
- if(!is.null(startIV1) ){ RTSRV1 = startIV1 }
- if(!is.null(startIV2) ){ RTSRV2 = startIV2 }
+ n = length(logprices1)
+ nbarK_cov = (n - K_cov + 1)/(K_cov)
+ nbarJ_cov = (n - J_cov + 1)/(J_cov)
+ adj_cov = n/((K_cov - J_cov) * nbarK_cov)
- if(is.null(startIV1)){ RTSRV1 = RTSRV( pdata1, noisevar = noisevar1 , K = K , J = J, eta = eta ) };
- if(is.null(startIV2)){ RTSRV2 = RTSRV( pdata2, noisevar = noisevar2 , K = K , J = J, eta = eta ) };
-
-
- I_K1 = 1*( logreturns_K1^2 <= eta*(RTSRV1 * vdelta_K + 2*noisevar1) );
- I_K2 = 1*( logreturns_K2^2 <= eta*(RTSRV2 * vdelta_K + 2*noisevar2) );
- I_J1 = 1*( logreturns_J1^2 <= eta*(RTSRV1 * vdelta_J + 2*noisevar1) );
- I_J2 = 1*( logreturns_J2^2 <= eta*(RTSRV2 * vdelta_J + 2*noisevar2) );
-
- if(eta ==9){ccc = 1.0415}else{ ccc = cfactor_RTSCV(eta=eta) }
-
- RTSCOV = adj * (
- ccc *(1/K)*sum(logreturns_K1*I_K1*logreturns_K2*I_K2)/mean(I_K1*I_K2)
- - ((nbarK/nbarJ) *ccc*(1/J)*sum(logreturns1*logreturns2*I_J1*I_J2 )/mean(I_J1*I_J2)));
- return(RTSCOV)
+ logreturns_K1 = logreturns_K2 = vdelta_K = c()
+ for (k in 1:K_cov) {
+ sel.avg = seq(k, n, K_cov)
+ logreturns_K1 = c(logreturns_K1, diff(logprices1[sel.avg]))
+ logreturns_K2 = c(logreturns_K2, diff(logprices2[sel.avg]))
+ vdelta_K = c(vdelta_K, diff(seconds[sel.avg])/secday)
+ }
+
+ logreturns_J1 = logreturns_J2 = vdelta_J = c()
+ for (j in 1:J_cov) {
+ sel.avg = seq(j, n, J_cov)
+ logreturns_J1 = c(logreturns_J1, diff(logprices1[sel.avg]))
+ logreturns_J2 = c(logreturns_J2, diff(logprices2[sel.avg]))
+ vdelta_J = c(vdelta_J, diff(seconds[sel.avg])/secday)
+ }
+
+ I_K1 = 1 * (logreturns_K1^2 <= eta * (RTSRV1 * vdelta_K + 2 * noisevar1))
+ I_K2 = 1 * (logreturns_K2^2 <= eta * (RTSRV2 * vdelta_K + 2 * noisevar2))
+ I_J1 = 1 * (logreturns_J1^2 <= eta * (RTSRV1 * vdelta_J + 2 * noisevar1))
+ I_J2 = 1 * (logreturns_J2^2 <= eta * (RTSRV2 * vdelta_J + 2 * noisevar2))
+ if (eta == 9) {
+ ccc = 1.0415
+ } else {
+ ccc = cfactor_RTSCV(eta = eta)
+ }
+ RTSCV = adj_cov * (ccc * (1/K_cov) * sum(logreturns_K1 * I_K1 *
+ logreturns_K2 * I_K2)/mean(I_K1 * I_K2) - ((nbarK_cov/nbarJ_cov) *
+ ccc * (1/J_cov) * sum(logreturns_J1 * logreturns_J2 * I_J1 *
+ I_J2)/mean(I_J1 * I_J2)))
+ return(RTSCV)
}
-TSCov = function (pdata,cor = FALSE, K = 300, J = 1, makePsd=FALSE)
+TSCov = function (pdata, cor = FALSE, K = 300, J = 1, K_cov = NULL, J_cov = NULL, makePsd = FALSE)
{
- #pdata is list with each item an xts object with 1 tick-by-tick timeseries
- # if(hasArg(data)){ rdata = data }
- # if(makeReturns){ rdata = makeReturns(rdata)}
-
- if(!is.list(pdata)){ n = 1 }else{
- n = length(pdata);
- if(n==1){pdata=pdata[[1]];}
- }
-
- if( n == 1 ){ return( TSRV( pdata, K = K , J = J ))}
-
- if( n > 1 ){
-
- cov = matrix(rep(0, n * n), ncol = n)
- diagonal = c()
- for (i in 1:n){
- diagonal[i] = TSRV(pdata[[i]], K = K , J = J )
- }
- diag(cov) = diagonal
- for (i in 2:n) {
- for (j in 1:(i - 1)) {
- cov[i, j] = cov[j, i] = TSCov_bi(pdata[[i]], pdata[[j]], K = K, J = J)
- }
- }
-
- if(cor==FALSE){
- if(makePsd==TRUE){cov = makePsd(cov);}
- return(cov)}
- if(cor==TRUE){
- invsdmatrix = try(solve(sqrt(diag(diag(cov)))),silent=F);
- if( inherits( invsdmatrix, "try-error") ){
- rcor = invsdmatrix%*%cov%*%invsdmatrix;
- if(makePsd==TRUE){rcor = makePsd(rcor);}
- return(rcor)}
- }
- }
+ if (!is.list(pdata)) {
+ n = 1
+ }
+ else {
+ n = length(pdata)
+ if (n == 1) {
+ pdata = pdata[[1]]
+ }
+ }
+ if (n == 1) {
+ return(RTAQ:::TSRV(pdata, K = K, J = J))
+ }
+ if (n > 1) {
+ cov = matrix(rep(0, n * n), ncol = n)
+ diagonal = c()
+ for (i in 1:n) {
+ diagonal[i] = RTAQ:::TSRV(pdata[[i]], K = K, J = J)
+ }
+ diag(cov) = diagonal
+ if( is.null(K_cov)){ K_cov = K }
+ if( is.null(J_cov)){ J_cov = J }
+ for (i in 2:n) {
+ for (j in 1:(i - 1)) {
+ cov[i, j] = cov[j, i] = RTAQ:::TSCov_bi(pdata[[i]],
+ pdata[[j]], K = K_cov, J = J_cov)
+ }
+ }
+ if (cor == FALSE) {
+ if (makePsd == TRUE) {
+ cov = makePsd(cov)
+ }
+ return(cov)
+ }
+ if (cor == TRUE) {
+ invsdmatrix = try(solve(sqrt(diag(diag(cov)))), silent = F)
+ if (!inherits(invsdmatrix, "try-error")) {
+ rcor = invsdmatrix %*% cov %*% invsdmatrix
+ if (makePsd == TRUE) {
+ rcor = makePsd(rcor)
+ }
+ return(rcor)
+ }
+ }
+ }
}
TSCov_bi = function (pdata1, pdata2, K = 300, J = 1)
More information about the Blotter-commits
mailing list