Consider the autoregressive model of order 2 $$X_{t}=\varphi_1X_{t-1}+\varphi_2X_{t-2}+\varepsilon_t,$$ where $\varepsilon_t$ are zero-mean normally distributed random variables with $\sigma^2$ variance such that these random variables are uncorrelated. Suppose that we have a sample for the above model with sample size, let say, somewhere between 100 and 1000.
As a part of a simulation I would like to know something about the probability distribution of the data. Some examples I made using MATLAB suggests that some normal distribution produces an excellent fitting, which - at least in my opinion - makes sense because the white noise process in the model is normally distributed.
However I am not skilled in the topic. Is there any result in the literature which can give me a theoretical base for this suggestion (that a sample of the above AR(2) model is normally distributed if considered as values of some random variable rather than as a time-series)? If there is not any, then how can I support the good fitting of some normal distribution?
At stationarity, the AR(2) model above is centered normal and its covariance structure $c_n=E(X_tX_{t-n})$ is given by $$ c_0=\frac{1-\varphi_2}{\Delta}\sigma^2,\qquad c_1=\frac{\varphi_1}{\Delta}\sigma^2, $$ where $$ \Delta=(1-\varphi_2)(1-\varphi_1^2-\varphi_2^2)-2\varphi_1^2\varphi_2, $$ and, for every $n\geqslant2$, $$ c_n=\varphi_1c_{n-1}+\varphi_2c_{n-2}. $$