Primero vamos a crear las variables de caminata aleatoria

set.seed(12345678)
y <- arima.sim(n=250, list(order=c(1,1,1), ma=0.25, ar=0.25), innov=rnorm(250))
x <- arima.sim(n=250, list(order=c(0,1,1), ma=0.65), innov=rnorm(250))

nota: el módelo ARIMA usa (p,d,q) para definir su orden. p= parte autorregresiva, d= Veces que se diferenció para q sea estacionaria, q= parte ma

Grafiquemos a ver como se ve:
layout(matrix(c(1, 2, 3), nrow = 3, ncol = 1, byrow = TRUE))
plot.ts(y, ylab=expression(y[t]), xlab="t", col="red", main = expression(y[t]))
plot.ts(x, ylab=expression(x[t]), xlab="t", col="red", main = expression(x[t]))
plot(x,y, ylab=expression(y[t]), xlab=expression(x[t]), col= "blue")

Usemos minímos cuadrados ordinarios

res1 <- lm(y~x)
r1 <- summary(res1)
attributes(r1)
## $names
##  [1] "call"          "terms"         "residuals"     "coefficients" 
##  [5] "aliased"       "sigma"         "df"            "r.squared"    
##  [9] "adj.r.squared" "fstatistic"    "cov.unscaled" 
## 
## $class
## [1] "summary.lm"
r1$r.squared
## [1] 0.7633916

Como podemos ver, parece ser que el 76% de la variación de y es explicada por x. Un r cuadrado muy alto. También la pendiente es significativa

Hagamos este proceso 10000 veces

N<- 10000
n<- 250
R2<-matrix(0,N,1)
for (i in 1:N) {
  y <- arima.sim(n=n, list(order=c(1,1,1), ma=0.25, ar=0.25), innov=rnorm(n))
x <- arima.sim(n=n, list(order=c(0,1,1), ma=0.65), innov=rnorm(n))
  R2[i]<-summary(lm(y~x))$r.squared
}
layout(matrix(c(1), nrow = 1, ncol = 1))
hist(R2,40, col = "green", main = "Histograma de los R cuadrados",
     ylab="Frecuencia", xlab="Valor R cuadrado")