Initialize AR(p) process by using Arima.sim

Hi I am currently trying to simulate an AR(4) process $y_t=0.67y_{t-1}-0.51y_{t-4}+\epsilon_t$ given that the initial value $y_1=1,y_2=2,y_3=3,y_4=4$ and $\epsilon_t\sim N(0,1)$. My code is given as below

> y <- arima.sim(34,model=list(ar=c(0.6,0,0,-0.5)),
start.innov=c(1,2,3,4),n.start=4,innov=c(0,0,0,0,rnorm(30,mean=0,sd=1)))

> y
Time Series:
Start = 1
End = 34
Frequency = 1
[1]  3.541600000  0.824960000 -1.785024000 -4.439014400 -5.947049868 -4.222318979 -2.326217678  1.144088356  5.150246659  4.159775404  5.767793661
[12]  0.967415814 -1.472489589 -3.311269208 -3.235976179 -1.787193697  0.855102209  1.619000292  3.488545520  2.288206839  1.247018920 -0.197260086
[23] -1.554881352 -2.202977785 -1.038750927 -1.255710320 -0.006073465  0.080668247  0.461419926 -0.274830999  0.806213750  0.216685532 -0.398132839
[34]  0.122593066


While what I am expecting is that the first four values of $y$ to be 1,2,3,4 say since I've set the first four entries of innov to be zero. Can someone help me with this problem? Thanks in advance.

• Welcome to CV! Hint: If you want to paste the block of code use the {} button. – Tim Dec 10 '14 at 7:46
• Do not crosspost between SO and CV: stackoverflow.com/questions/27396227/… , posting your question on one site is enough. – Tim Dec 10 '14 at 8:38

You have the following system:

\begin{eqnarray} \left\{ \begin{array}{lcl} y_5 &=& 0.67 y_4 - 0.51 y_1 + \epsilon_5 \\ y_6 &=& 0.67 y_5 - 0.51 y_2 + \epsilon_6 \\ y_7 &=& 0.67 y_6 - 0.51 y_3 + \epsilon_7 \\ y_8 &=& 0.67 y_7 - 0.51 y_4 + \epsilon_8 \end{array} \right. \end{eqnarray}

with $\epsilon_t \sim NID(0, 1)$. In order to achieve $y_5=1,\, y_6=2,\, y_7=3\,$ and $y_8=4$, you can set the first eight innovations $\epsilon_t$ to zero and solve the system above for $y_1, y_2, y_3$ and $y_4$ given the desired values for $y_5, y_6, y_7$ and $y_8$. That is, the system to solve becomes, in matrix form:

\begin{eqnarray} \left( \begin{array}{cccc} -0.51 & 0 & 0 & 0.67 \\ 0 & -0.51 & 0 & 0 \\ 0 & 0 & -0.51 & 0 \\ 0 & 0 & 0 & -0.51 \end{array} \right) \left( \begin{array}{c} y_1 \\ y_2 \\ y_3 \\ y_4 \end{array} \right) = \left( \begin{array}{c} 1 \\ 1.33 \\ 1.66 \\ 1.99 \end{array} \right) \,. \end{eqnarray}

This gives the solution: $y_1=-7.086890,\, y_2=-2.607843,\, y_3=-3.254902\;\,$ and $\,y_4=-3.901961$.

The vector $(1, 1.33, 1.66, 1.99)$ is obtained as follows: the first element is $y_5$, for which we want the value $1$ ($\epsilon_5$ is set to zero); the second element is $-0.51y_2 = y_6 - 0.67y_5 = 2 - 0.67\times1 = 1.33$ (the desired value for $y_6$ is $2$ and $\epsilon_6$ is set to zero); the third element is $-0.51y_3 = y_7 - 0.67y_6 = 3-0.67\times 2 = 1.66$; and from the last equation $-0.51y_4 = y_8 - 0.67y_7 = 4 - 0.67\times 3 = 1.99$.

Now, upon this result, you should define the arguments n.start, start.innovand innov that are passed to arima.sim. A similar example is given here. For this case I couldn't figure out the right definition of these arguments, there may be something else done by arima.sim that I am overlooking. Nevertheless, you can still generate your data as follows:

set.seed(123)
x <- rep(0, 20)
x[1:4] <- c(-7.086890, -2.607843, -3.254902, -3.901961)
eps <- c(rep(0, 8), rnorm(12))
for (i in seq(5,20))
x[i] <- 0.67 * x[i-1] - 0.51 * x[i-4] + eps[i]
x[-seq(4)]
#  [1]  1.000  2.000  3.000  4.000  1.610 -0.172 -0.086 -2.027 -2.050  0.429
# [11]  0.793  0.300  0.560 -0.290  0.626  0.626


After the auxiliary observations $y_1$ to $y_4$ that were found above, the series continues with the desired values $1, 2, 3, 4$.