How to understand the square of an AR(1) process?

Question

I generate an AR(1) process as follows:

x=arima.sim(list(order = c(1,0,0),ar=0.67),n=1000,sd=sqrt(0.55))

When I square it, and fit AR(1) to the squared process, it seems still to be fitted well using AR(1) model:

y=x^2
fit=ar(y,order.max=1, method="ols")
acf(fit$resid[2:1000])
Box.test(fit$resid,lag=round(log(length(y))),type="Ljung-Box",fitdf=1)$p.value

The ACF plot of the model (for the squared process) residuals and the large residual test p value (larger than .05) both show that the residuals are independently distributed. What does this indicate?

One simulation of time series $x$ and $y$ are as follows:

Answer 1

I like the question. I think I have convinced myself of why you are seeing this, but it may be wishful thinking.

The acf of the residuals shows that the residuals have low autocorrelations at lags 1 and above. It does not mean that they are independent. I think I can see why these autocorrelations are likely to be small. You can write the original process as $$x_i = \alpha x_i + \varepsilon_i$$ and squaring gives $$x_i^2 = \alpha^2 x_i + 2\alpha \varepsilon_i x_i + \varepsilon_i^2$$ Assuming that the ar() function gives the right value of $\alpha^2$, the residuals are $$r_i = 2\alpha\varepsilon_i x_i + \varepsilon_i^2$$ but $x_i = \sum_{k=0}^\infty \alpha^k \varepsilon_{i-k}$ because this is an AR(1) process, so $$r_i = 2\alpha \varepsilon_i^2 + 2\alpha^2 \varepsilon_{i-1} \varepsilon_i + 2\alpha^3 \varepsilon_{i-2}\varepsilon_i + \cdots + \varepsilon_i^2$$ Now look at lag 1 (it's similar for higher lags). If we try to calculate the covariance of $r_i$ with $r_{i-1}$ and use the linearity of covariance, we end up with possibly nonzero terms of the form $$cov(\varepsilon_{i-1}^2, \varepsilon_i \varepsilon_{i-1})$$ and $$cov(\varepsilon_{i-1} \varepsilon_j, \varepsilon_{i} \varepsilon_j)$$ with other terms dropping out due to independence. But both of these are likely to be small (as can be checked by simulation) because if you imagine plotting $\varepsilon_{i-1}^2$ versus $\varepsilon_{i-1}$ you will get a parabola, and multiplying the points by an independent random amount in the y-direction will give you something which is likely to have a roughly horizontal line of best fit, hence small correlation, hence the covariance should be close to zero. Similarly, plotting $\varepsilon_{j}$ versus $\varepsilon_{j}$ and then multiplying all the points by independent random amounts in the x and y-directions will give something that looks very much like a random scatter of points, so the second covariance should be small. Again, it can be checked via simulation.

Notice that this explanation indicates that you should see the same behaviour regardless of the value of sd, which does indeed seem to be the case (I tried it with sd=1 and sd=5).

I do not think the same explanation works for exponential white noise. I tried doing it with exponential white noise and sometimes the software fits and AR(1), sometimes an AR(4).

Answer 2

Quick and dirty - my thoughts on this. This might give you something to go on though I dunno if it is helpful. All usual assumptions

$$x_i=\rho x_{i-1} + \epsilon_t \\ y_i=x^2\\ y_i=(\rho x_{i-1} + \epsilon_t)^2\\ y_i=\rho^2 x_{i-1}^2 + 2\rho x_{i-1} \epsilon_t + \epsilon_t ^ 2\\ y_i=\rho^2 x_{i-1}^2 + 2\rho x_{i-1} \epsilon_t + \delta\\ E[\delta] = E[\epsilon_t^2] = \sigma^2 $$

Okay so then you get for some estimator which uncorrectly assuming AR(1) without including the quadratic thing

$$y_i^*=\phi^*y_{i-1}\\ \delta^* = y_i-y_i^*\\ E[\delta^*]=E[ y_i-y_i^*]=E[\phi x_{i-1}^2 + 2\rho x_{i-1} \epsilon_t + \delta - \phi^*y_{i-1}]\\ = (\phi-E[\phi^*]) y_{i-1} + \sigma^2 $$

So assuming the estimator is somehow unbiased we have $E[\delta^*] = \sigma^2$

$$Cov[\delta^*_t, \delta^*_{t-1}] = E[(\delta^*_t-E[\delta^*_t])(\delta^*_{t-1}-E[\delta^*_{t-1}])]\\ = E[(\delta^*_t-\sigma^2)(\delta^*_{t-1}-\sigma^2)] \\ = E[\delta^*_t \delta^*_{t-1} - \sigma^2 \delta^*_{t-1} - \sigma^2 \delta^*_{t} + \sigma^4]\\ = E[\delta^*_t \delta^*_{t-1}] = \sigma^4$$

I can see where you are coming from, but maybe this is just too small a number? 0,55^4 is around 0,09. What is the result for a much larger simulated standard deviation?

Maybe this also breaks down because the estimator is no longer unbiased. I will continue this for the biased case later but maybe you can go with this and check my thoughts there's probably some type of error in there.