I have created a VAR model in R (using the command VAR) and my model reports an R-squared of about 0.8 for the variable I'm most interested in. I'm trying to replicate that result by calculating the R-squared by hand. I have extracted the fitted values (using the command fitted(model_name)) and then I have used the following equation:

$$R^2 = 1- { \sum_{} (y-\hat y)^2 \over \sum_{}(y - \overline y)^2}.$$

With $y$ = actual value, $ \hat y$ = predicted (fitted) value and $ \overline y$ = the average of the actual value.

I am taking into account the fact that my model uses two lags, so I have two fitted values less than the total amount of points used to train. I'm making sure to pair the actual and predicted values correctly.

I have tried using this equation before in regular linear regression and I am able to replicate the value I get from R, however I cannot do it with this VAR model. Is it because when the R-squared is calculated for time series, something different is done? I have tried googling this and cannot find any reference to a different equation for the case of time series.

Any comments or suggestions as to what could be wrong will be greatly appreciated.

  • $\begingroup$ maybe the mean is also only calculated from the "clipped" time series after discarding initial observations due to the lags? $\endgroup$ – Christoph Hanck Jul 7 '20 at 7:59
  • $\begingroup$ That is a great observation. I tried that as well and it did not work $\endgroup$ – amilkar0417 Jul 7 '20 at 19:35
  • $\begingroup$ maybe you can post some code of what you did - please see my answer below, which suggests that it does work $\endgroup$ – Christoph Hanck Jul 8 '20 at 5:45

In principle, the procedure for computing $R^2$ should indeed be as you describe also for VAR models, as these are also simply fitted by OLS.

The discrepancy indeed seems to be in how VAR computes the numerator in the formula for $R^2$.

Consider the code below:


# generate some data
n <- 10
y1 <- rnorm(n)
y2 <- rnorm(n)
y <- cbind(y1,y2)

# R^2 from vars package
var.model <- VAR(y, p = 1)
R2.vars.1 <- summary(var.model$varresult$y1)$r.squared

# handmade
u1 <- resid(var.model)[,1]
# gives different results when computing numerator from full series:
#R2.1 <- 1-sum(u1^2)/sum((y1-mean(y1))^2) 

# numerator of R^2 when only using the portion that can actually serve as a dependent variable
R2.1 <- 1-sum(u1^2)/sum((y1[2:n]-mean(y1[2:n]))^2) 

all.equal(R2.1,R2.vars.1) # TRUE (identical and == say FALSE because of small numerical differences)

Your Answer

