Calculating the means in a for loop in R-CodePudding

Good afternoon all,

I am trying to find the standard prediction error of a time-series that I generate through a simulation run which is defined through the function called sim_11 with 250 simulations. This is provided in the first batch of code below.

The second batch creates a time-series model (AR(1)) and tries to predict the next 5 values, and I do this overall 250 simulations. For each simulation, I should be able to get 5 prediction errors and after 250 simulations I should have a resulting table of 250 rows and 5 columns. However, when I try to set this up in the for loop, I end up with only 250 single values when in fact I should end up with a 250 by 5 table/matrix. I believe the error to be in the

pred_error_AR1_100[i]<-table((pre_AR1_100$se[1]),(pre_AR1_100$se[2]),
                           (pre_AR1_100$se[3]),(pre_AR1_100$se[4]),
                           (pre_AR1_100$se[5]), ncol=5)

part however I am not able to figure out where or what the format should be.

Thank you in advance.

The two code batches are provided below for replication.

# Setup the simulation run with 100 observations and 250 simulations
sim_11=function(){
  e<-rnorm(200, mean=0, sd=0.2) # Produces 200 white noise values
  Y_t=c(0,0)  # Fills in the first 2 observations as a lag of 2 can be handled
  for (i in 3:length(e)){
    f1<- 0.138 (0.316 0.982*Y_t[i-1])*exp(-3.89*(Y_t[i-1])^2)
    f2<- -0.437-(0.659 1.260*Y_t[i-1])*exp(-3.89*(Y_t[i-1])^2)
    Y_t[i]<-f1*Y_t[i-1] f2*Y_t[i-2] e[i]
  }
  Y_t<-Y_t[101:200] # Removes the first 100 observations
  Y_t # Prints the 100 observations
}

lapply(1:250, function(x) sim_11()) # Provides the results of the 250 simulations
x_100_lstar=replicate(250,sim_11()) # Places all results into one matrix

pred_error_AR1_100=0
# controls<-list(gammaInt=c(0.1,2000), nGamma=50)
for (i in 1:ncol(x_100_lstar)){
  AR1_100<-ar(x_100_lstar[,i])
  pre_AR1_100<-predict(AR1_100, n.ahead=5)
  pred_error_AR1_100[i]<-table((pre_AR1_100$se[1]),(pre_AR1_100$se[2]),
                           (pre_AR1_100$se[3]),(pre_AR1_100$se[4]),
                           (pre_AR1_100$se[5]), ncol=5)
}
pred_error_AR1_100

CodePudding user response：

To get your loop to work, you need to initialize pred_error_AR1_100 as an n-by-5 matrix, then modify the rows one at a time. You should not be using table here. See ?matrix and ?Extract for details on constructing, accessing, and modifying matrices.

n <- ncol(x_100_lstar)
pred_error_AR1_100 <- matrix(NA, n, 5)
for (i in seq_len(n)) {
  AR1_100 <- ar(x_100_lstar[, i])
  pre_AR1_100 <- predict(AR1_100, n.ahead = 5)
  pred_error_AR1_100[i, ] <- pre_AR1_100$se
}

In these situations, though, it is safer and faster to use apply than to write the loop yourself:

## Here, 'x' represents the result of one realization of 'sim_11()'
f <- function(x) {
  AR1_100 <- ar(x)
  pre_AR1_100 <- predict(AR1_100, n.ahead = 5)
  pre_AR1_100$se
}

## Apply function 'f' to each column of 'x_100_lstar'
pred_error_AR1_100 <- t(apply(x_100_lstar, 2, f))

In the last line, the result of apply, a 5-by-n matrix, is transposed to get an n-by-5 matrix.

FWIW, sim_11() would be slightly faster if you initialized Y_t as a vector of length 200, like so:

Y_t <- rep.int(NA, 200)
Y_t[1:2] <- 0

rather than incrementing the length by 1 in each iteration.