Problem: When trying to calculate the variance of timeseries sums I get a negative variance, mostly due to autocovariances at large lag steps. Does not seem realistic.
I have a timeseries which is calculated from another timeseries using a regression equation. I would like to propagate the uncertainty in the regression to the final timeseries. Then I want to sum (or take mean values) different segments of the timeseries over different timeperiods, and get the uncertainty of the sums. The timeseries is originally in 1 hour frequency and I want to sum over periods of 1 day (resampling to daily frequency) up to several years. The timeseries is strongly autocorrelated at short lag times.
For getting the variance of the sum (in the case of 3 elements being summed): $$Var(a+b+c)= \\ Var(a)+Var(b)+Var(c) + 2 \times (Cov(a,b) + Cov(a,c)+Cov(b,c))$$
I use r
for the calculations. I get the variances for each timeseries element as $SE^2$, where $SE$ is the standard error (se.fit
) returned from r's predict()
function using the regression model. The covariances I get from the autocovariance function acf()
.
Here is some code and a selection of the data (excuse clumsy R code, I'm very new to R):
#tsY is the predicted timeseries from the regression
tsY=c(81.4, 79.0, 83.4, 81.7, 75.7, 68.3, 62.3, 57.2, 52.6, 48.8, 45.4, 42.6, 39.9, 37.6, 35.6, 33.8, 32.2, 30.8, 29.6, 28.4, 27.3, 26.2, 25.0, 23.9)
#tsSE is the standard error from the prediction (se.fit)
tsSE=c(1.55, 1.49, 1.60, 1.56, 1.41, 1.23, 1.09, 0.97, 0.87, 0.78, 0.71, 0.65, 0.60, 0.55, 0.51, 0.48, 0.45, 0.42, 0.40, 0.38, 0.36, 0.34, 0.32, 0.30)
tsVar=tsSE^2
#create a matrix of the autocovariances at different lag times, diagonal is lag=0
#rows and columns are indicies in timeseries
covmat<-matrix(numeric(0), length(tsY),length(tsY))
for ( i in (1:(length(tsY)) ) ) {
if (i == 1) {
autocov<-acf(tsY, type='covariance', lag.max= length(tsY))
autocovvec<-autocov$acf[1:nrow(autocov$acf)]
covmat[i:length(tsY),i]=autocovvec
} else {
autocov<-acf(tsY[-(1:i-1)], type='covariance', lag.max= length(tsY))
autocovvec<-autocov$acf[1:nrow(autocov$acf)]
covmat[i:length(tsY),i]=autocovvec
}
}
# sum the matrix columns, but not the diagonal
sumofColumns <- rep(NA, ncol(covmat))
for (i in (1:ncol(covmat))) {
if (i == 1) {
sumofColumns[i]=sum(covmat[-(1),i])
} else{
sumofColumns[i]=sum(covmat[-(1:i),i])
}
}
sumofCov=sum(sumofColumns) # sum of the covariance (Cov(a,b) + Cov(a,c)+...)
sumofVar=sum(tsVar) # sum of the variances of each timeseries element
varofSum=sumofVar+2*sumofCov # variance of the sum of the timeseries
# from the covmat the negative variance occurs at larger lag times.
acf(tsY, type='covariance', lag.max= length(tsY))
> sumofCov
[1] -1151.529
> varofSum
[1] -2283.246
So I have the following questions:
Did I completely misunderstand how to calculate variance of sums?
Is it better to use a cutoff from the max lags to be considered in the autocovariance? If so how would one determine this? This would especially be important with the complete data where the length is several thousand.
3. Why is the covariance negative in this sample data at large? When plotting tsYBecause it is the variation in direction from their means.plot(tsY)
it looks like the covariance/correlation should remain positive.
EDIT:
Comment on question 2 above: I have realized that using n-1 lags, as above in the code, does not make a lot of sense. There appear to be few different ways to determine the maximum lags to consider. Box & Jenkins (1970) suggest n/4 and R by default 10*log10(n). This does not answer the question however, of how to determine an appropriate cutoff for summing the covariances.
Does it make sense to look at the partial autocorrelation (function pacf()), in order not to overestimate the effect of the auto covariance in the summation term? The partial autocorrelation for my data is significantly different from zero only at 1 or 2 lags. Similarly, fitting an AR model using ar() function, I also get an order of 1 or 2.
Cheers
Related post Variance on the sum of predicted values from a mixed effect model on a timeseries