"Any* accuracy measure is likely to be worse out-of-sample than in-sample, because we fit the model *to the training sample*. If we fit it too closely (fitting noise), this results in so-called overfitting.

However, even if we don't overfit and get the data generating process (DGP) just right, it is still extremely likely that the DGP evolves slowly over time, so the out-of-sample DGP is different from the one we just assumed we got right in fitting. Result: worse out-of-sample than in-sample performance.

Incidentally, [$R^2$ is not commonly used in assessing forecast accuracy](https://stats.stackexchange.com/q/195387/1352), even though [it is equivalent to MSE](https://stats.stackexchange.com/a/327482/1352) (which *is* a common accuracy measure).

EDIT: [you ask whether this is more of a coincidence or of a rule](https://stats.stackexchange.com/questions/427822/out-of-sample-and-in-sample-forecasting-r-squared/427823?noredirect=1#comment798145_427823). 

Well, we do fit to historical observations, not to out-of-sample ones. And I think it's fairly uncontroversial that "DGP drift" exists. So I'd say it's a rule. So much so that I would be very surprised at any forecasting method that systematically *didn't* exhibit it. (Of course, single counterexamples will always exist, simply because of random variability.)

As an example, let's run through the first 2,100 time series of the M3 dataset. In each case, we fit a model using `auto.arima()` and check whether the in-sample $R^2$ is smaller than the out-of-sample $R^2$ (which, per above, is equivalent to the MSE, so by monotonicity, to the RMSE).

The result: the in-sample $R^2$ is less than or equal to the out-of-sample one in 1,536 series (and the other way around in 564 series). Which is actually fewer than I had expected.

R code:

    library(forecast)
    library(Mcomp)
    
    in_sample_leq_out_of_sample <- rep(NA,length(M3))
    
    pb <- winProgressBar(max=length(M3))
    for ( ii in seq_along(M3) ) {
    	setWinProgressBar(pb,ii,paste(ii,"of",length(M3)))
    	foo <- accuracy(forecast(auto.arima(M3[[ii]]$x),h=length(M3[[ii]]$xx)),M3[[ii]]$xx)
    	in_sample_leq_out_of_sample[ii] <- foo[1,2]<=foo[2,2]
    }
    close(pb)
    
    table(in_sample_leq_out_of_sample)