Nested Loop for summation

[FactOREO]

Hello community,

I've started working with R just a few weeks ago in university and was wondering about nested loop functions with multiple sums. In particular I have an ARMA(2,4) process, which I want to simulate. I found all neccesary parameters and first simulated the time series with a simple for loop:
for (i in 5:length(t)) {
Y[i] = a + phi.1 * Y[i-1] + phi.2 * Y[i-2] + Z[i] + theta.4 * Z[i-4] + theta.3 * Z[i-3] +theta.2 * Z[i-2] + theta.1 * Z[i-1]
}
where phi.1 and so on are the parameters I defined before. Because of the lag up to [i-4] I start the loop at 5. t is the time vector with 1,000 values and Z[i] is just white noise.
Since this is a pretty sh.. code, I tried to simplify the summation for phi and theta with 2 other for-loops, but unfortunately it doesn't work in the slightest. I think I do understand why it is not working the way Iam trying it, but I wasn't able to come up with a suitable soluation until now.
My try:
for (i in 5:length(t)) {
for (j in 1:2){
for (k in 1:4){
Y[i] = sum(phi[j]*Y[i-j]) + Z[i] + sum(theta[k]*Z[i-k])
}
}
}
I think it is not equivalent because R tries to do the k-loop (1:4) at i=5 and j=1, then goes to i=5 and j=2 and repeats k=1:4. But I want j=1:2 and k=1:4 at the same time for i=5.

Can anyone give me at least a hint, how I can handle this procedure, so that R computes the AR and MA part simultaneously for every i-th value I try to compute?

Sincerely

FactOREO

[technocrat]

See the FAQ: How to do a minimal reproducible example reprex for beginners. The idea is to have a cut and paste example that will load all needed libraries and data to illustrate the problem. Without it, only general guidance can be offered.

The following example shows stepwise development of the solution to this type of problem. Below, it was easy to see that with the presence of y[x-2] in the function that the for loop would have to start with 3. I suspect a similar approach will allow you to get the logic right in your case.

# Simulate y
y <- numeric(100)
e <- rnorm(100)
for(i in 2:100)
  y[i] <- 0.6*y[i-1] + e[i]

# define constants

a <- pi
phi.1 <- 0.6
phi.2 <- 0.2

# hardwired function to make illustration easier

make_term <- function(x) a + phi.1 * y[x-1] + phi.2 * y[x-2] 

# illustrate application

holder <- list()
for (i in 3:length(y)) holder[i] = make_term(i)
unlist(holder)
#>  [1] 2.6234268 2.3265332 2.0950606 2.7646958 2.8868463 2.9229873 2.9721795
#>  [8] 4.1038078 4.9533835 4.8940731 4.7022580 4.3387726 3.9995287 2.7366209
#> [15] 2.4299101 1.5431100 0.3188871 0.8187863 2.1499154 2.6044797 2.1201821
#> [22] 1.4736645 2.7649998 3.1557277 2.7981470 2.1823276 2.5958512 3.1551375
#> [29] 4.0576772 3.7542651 3.5927147 3.0292660 3.7635644 3.3224826 3.1727496
#> [36] 3.5903654 3.6305808 3.7284071 2.3771208 2.1907763 1.9172209 1.4377688
#> [43] 2.3065391 3.0730440 3.3824714 3.5517975 3.0713545 3.9001606 3.5860109
#> [50] 3.4592384 3.0830691 3.1155525 3.2328833 3.6112411 3.7630496 2.5719653
#> [57] 2.6623233 2.6805534 3.2577633 3.1121410 3.0978636 3.4197126 3.8561072
#> [64] 3.4222195 4.0599125 5.1739849 4.3072138 3.4602369 2.9775378 2.2128272
#> [71] 1.9474705 2.7053496 3.3890641 3.6145090 3.4878419 2.5550897 3.4021988
#> [78] 3.4250704 4.3320918 3.7166646 3.3404745 3.3253659 3.8951965 4.6573359
#> [85] 4.6854418 3.7268465 2.7618039 2.2864325 3.1298235 4.5472345 4.9670926
#> [92] 4.3350595 3.8169973 3.6981643 3.6270248 3.6055920 3.4311680 2.7240743

[FactOREO]

Thanks for your comment. Iam going to provide the reprex right down here. Sorry for my mistake.

t <- seq(0.01, 10, 0.01)
Z = rnorm(length(t), 0, 1)

#coefficients
#intercept
a <- -0.005314   

#AR-coefficients
phi.1 <-  0.95038       
phi.2 <- -0.25543         
#MA-coefficients
theta.1 <- 0.04334          
theta.2 <- 0.08012       
theta.3 <- -0.07681     
theta.4 <- -0.01477 

#ts construction
Y <- vector()
Y[1:4] = Z[1:4]
for (i in 5:length(t)) {
Y[i] = a + phi.1 * Y[i-1] + phi.2 * Y[i-2] + Z[i] + theta.4 * Z[i-4] + theta.3 * Z[i-3] +theta.2 * Z[i-2] + theta.1 * Z[i-1]
 }
 acf(Y)

pacf(Y)

The given code reproduces a time series with acf cutting at lag=4 and pacf cutting at lag=2.

And there is no way to get arround using a function resp. a list to summarize the sums? I was hoping to just use a nested loop to be able to reduce the code and also summarize the variables in the form of

phi <- c( 0.95038  , -0.25543 )
theta <- c(0.04334, 0.08012, -0.07681, -0.01477)

and then adress the components of the parameter-vectors one by the other.

Thanks for your help and kind regards

FactOREO

[technocrat]

@FactOREO

Sorry for my mistake.

Not at all. Good to have the starting point. Whenever I get a bit of picky code like that I like to put it in a function, if for no other reason than to be able to get back to a safe place

remake_Y <- functxon(x) {
  Y[x] = a + phx.1 * Y[x-1] + phx.2 * Y[x-2] + Z[x] + theta.4 * Z[x-4] + theta.3 * Z[x-3] +theta.2 * Z[x-2] + theta.1 * Z[x-1]
}
### snip
for (i in 5:length(t)) remake_Y(i)

And I would do the same for the triple loop's function

sum_phi <- function(i,j) sum(phi[j]*Y[i-j]) 

sum_theta <- function(i,k) sum(theta[k]*Z[i-k])

make_again <- function(i,j,k) Y[i] = sum_phi(i,j) + Z[i] + sum_theta(i,k)

The purpose is to keep shortening the for loop to find what happens at a particular place. Now, since I'm approaching this not as an ARMA(2,2) analysis but as an abstract question, I now have to ask, what is the relationship between \phi in the last function and \phi_1 and \phi_2 in the first?

[FactOREO]

Am I right that you are asking, what the relationship between

@technocrat

ϕ1\phi_1 and ϕ2

-> phi.1 and phi.2 in my code
and

@technocrat

ϕ

-> phi <- c(phi.1, phi.2) in my code
is? If that's the case, than it's just the combination of phi 1 and 2 for the autoregressive part of the ARMA time series in one vector.

Otherwise you may ask the question again in a different form, since I am not sure if I understood it correctly

[technocrat]

Questions are the hardest part of data science.

In the first loop we were separately dealing with the two phi constants as scalars. Do you mean in the second loop that the two phi constants will be combined as a vector? ("Scalar" as in an R vector of length one.)

[FactOREO]

Yes, exactly. In my second try (the triple loop) I tried to combine the two scalars (e.g. phi1 and 2) in one 2x1 vector with components phi1 and 2. The loop should then pick for every i-th value of the time series Y(t) both phi values (as well as theta-value from 1 to 4) and match them with the corresponding data from the past (e.g. Y(t-1) and Y(t-2) for phi 1 and 2).
Therefore there is no other connection between the scalars phi 1 and 2 and the vector phi, than that phi consists of phi 1 and 2. I hope, this answered your question appropriately.