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

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:

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: