How to make the previous loop value effective in current iteration in "apply" ("apply"-family version of a "for" loop)?

Question

for loop:

N <- 500; 
ro <- 0.6; a <- 1
set.seed(1)
v <- ts(rnorm(N,0,1))

# [1] -0.626453811  0.183643324 -0.835628612  1.595280802  0.329507772 
# [496] -1.108909998  0.307566624 -1.106894472  0.347653649 -0.873264535

y <- ts(rep(0,N)) # y[1]:=0 defined
for (t in 2:500){ y[t] <- a + ro*y[t-1] + v[t] }

y
# [1]  0.00000000  [2] 1.18364332  
# [3] 0.87455738=1+0.6*1.18364332+(-0.835628612)  
# [4] 3.12001523=1+0.6*0.87455738+(1.595280802)
# [499]  2.55513301  1.65981527

mean(y) #2.549763

I wanted to convert the above for loop to an apply-family version:

N <- 500; ro <- 0.6; a <- 1
set.seed(1)
v <- as.data.frame(rnorm(N,0,1)) 
# 1: -0.626453811 2: 0.183643324 3: -0.835628612 4: 1.595280802 5: 0.329507772
# 496: -1.108909998 497: 0.307566624 498: -1.106894472 499: 0.347653649 500:-0.873264535
y <- as.data.frame(c(y1=0,rep(0,N-1))) # index starts at 1.  y[1]:=0 defined
y <- c(y[1,], unlist(apply(as.matrix(2:500), 1, function(t) { y[t,] <- a + ro*y[t-1,] + v[t,] })))
y
# [1]  0.000000000  [2] 1.183643324  
# [3] 0.164371388=1+0.6*0+(-0.8356286124) 
# [4] 2.595280802=1+0.6*0+(1.5952808025)
# [496] -0.108909998  1.307566624 -0.106894472  1.347653649  0.126735465

mean(y) # 1.021897

The above code does not give the results in for loop. I discovered what is wrong: in apply, in the iteration equation, previous values of y are not used; instead, y <- as.data.frame(c(y1=0,rep(0,N-1))); i.e. yt=0 is used for all t.

What to do to make successful apply-family?

`y[-1] <- ro^(2:N) * cumsum((a+y[1]+v[-1]) / ro^(2:N))` may avoid an explicit `for` loop — Henry, Mar 23 '18 at 09:40
@Cath, `apply` seems to be smarter to me. That said, I just learnt `apply` is _not_ a vectorized solution alternative for `for` loop. — Erdogan CEVHER, Mar 23 '18 at 11:20
@Henry, Thx a lot for your very elegant solution. It seems to me that your solution is stemmed from Wold representation. I will write it as one alternative solution to the problem (of neating `for` solution). I will add equaivalent math latex expression when I figured out. — Erdogan CEVHER, Mar 23 '18 at 11:39
@Henry, could you please write the math equivalent of `y[-1] <- ro^(2:N) * cumsum((a+y[1]+v[-1]) / ro^(2:N))` or give a reference that leads you to write that one-liner? I could not figure out the math equivalent. — Erdogan CEVHER, Mar 24 '18 at 18:40
@ErdoganCEVHER: the cumulative sum is based on `y[m]/ro^m - y[m-1] = a/ro^m +v[m]/ro^m` and needing to incorporate `y[1]` which on review I have may misplaced - perhaps something like this revision `y[-1] <- ro^(2:N) * (cumsum((a+v[-1]) / ro^(2:N)) + y[1]/ro)` could do that better — Henry, Mar 24 '18 at 19:56

score 2 · Answer 1 · edited Mar 23 '18 at 11:59

2

Well here is a recursive function which does the same as your loop:

my.fct <- function(t, initial.value){
  if(t == 1){
    return(initial.value)
  }
  if(t > 1){
    vec <- my.fct(t-1, initial.value)
    val <- vec[length(vec)]
    newval <- a + ro*val + v[t]
    newvec <- c(vec, newval)
    return(newvec)
  }
}
# Indeed it yields the same result:
> sum(my.fct(N, 0) == y)-N
[1] 0

And you can always pass a; ro; v as arguments as well:

N <- 500; ro <- 0.6; a <- 1
set.seed(1)
v <- ts(rnorm(N,0,1))
my.fct(500,0) # yields exactly the same values in the "for" loop

Comparing performances:

# function
start.time <- Sys.time()
y <- ts(my.fct(N, 0))
end.time <- Sys.time()
time.taken <- end.time - start.time
time.taken
Time difference of 0.006044865 secs

# loop
start.time <- Sys.time()
y <- ts(rep(0,N))
for (t in 2:500){ y[t] <- a + ro*y[t-1] + v[t] }
end.time <- Sys.time()
time.taken <- end.time - start.time
time.taken
Time difference of 0.01205611 secs

edited Mar 23 '18 at 11:59

Erdogan CEVHER

1,788
1
21
40

answered Mar 23 '18 at 09:18

niko

5,253
1
12
32

2

Yes in that case you should increase the max number of nested expressions `options(expressions=...)` and in some instances as well the size of the protection stack (if the error `protection stack overflow` appears). See here for more details https://stackoverflow.com/questions/22003021/explanation-of-r-optionsexpressions-to-non-computer-scientists – niko Mar 23 '18 at 09:48
2

@Cath R has an option to limit the number of nested expressions (as suggested `options(expresssions=)`: 'expressions': sets a limit on the number of nested expressions that will be evaluated. Valid values are 25...500000 with default 5000... – Eumenedies Mar 23 '18 at 09:48
`my.fct(6000,0)` resulted in `Error: evaluation nested too deeply: infinite recursion / options(expressions=)?` with my default expressions=5000. It reduced my eagerness to accept the above solution as an answer. – Erdogan CEVHER Mar 23 '18 at 12:10
`options()$expressions` should be at least as high as argument `t`, so try setting `options(expressions = 50000)`. Also, as mentioned, if `protection stack overflow` appears then you'll have to increase `--max-ppsize` (c.f. https://stackoverflow.com/questions/28728774/how-to-set-max-ppsize-in-r) – niko Mar 23 '18 at 12:16

Erdogan CEVHER · Accepted Answer · 2018-04-17T07:58:39.447

AR(1) process:

$$ y_t= \alpha + \rho y_{t-1} + \nu_t, \hspace{1cm} t=2,...,n+1 $$ $$ (\nu_t \sim N(0,1)) $$

(t can be finalized at t=n if one includes initial $$ y_1 $$ value in N obs)

By math'l equivalent: $$ y_t=ρ^t (\frac{y_1}{\rho} + \sum_{j=2}^{t} \frac{\alpha+\nu_{j}}{\rho^j} ) (t=2,…,n; y_1 is prespecified) $$

Solution 1 (Henry's elegance):

rm(list = ls()) # Clear workspace by deleting all objects
N <- 500; ro <- 0.6; a <- 1
set.seed(1)
v <- ts(rnorm(N,0,1))
y <- ts(rep(0,N)) # y[1]:=0 defined
y[-1] <- ro^(2:N) * (y[1]/ro + cumsum((a+v[-1]) / ro^(2:N)))
y # yields exactly the same values in the "for" loop 
(mean(y)) # 2.549763

$$\rho^t\rightarrow 0$$ as $$ t \rightarrow \infty $$, even very early: 0.6^1458=4.940656e-324 ; 0.6^1459=0 This blows the model. When N=1400, it gives y[1388]=1.114; 1389...1396:Inf; 1397...1400:NaN

Solution 2 (Nate's formalism):

rm(list = ls()) #  Clear workspace by deleting all objects
my.fct <- function(t, initial.value){
  if(t == 1){ return(initial.value) }
  if(t > 1){
    vec <- my.fct(t-1, initial.value)
    val <- vec[length(vec)]
    newval <- a + ro*val + v[t]
    newvec <- c(vec, newval)
    return(newvec)
  }
}

N <- 500; ro <- 0.6; a <- 1
set.seed(1)
v <- ts(rnorm(N,0,1))
my.fct(500,0) # yields exactly the same values in the "for" loop
mean(my.fct(500,0)) # 2.549763

N=833 works. N=834 gives Error: evaluation nested too deeply: infinite recursion / options(expressions=)? with the default options(expressions=5000).

(I could not figure out how to make math and text in-line in latex above)

How to make the previous loop value effective in current iteration in "apply" ("apply"-family version of a "for" loop)?

2 Answers2