Eliminate the need for a loop in a simple example-CodePudding

I have a problem that I am currently solving using a loop, but something tells me that it is possible to do so without. The main reason that I think it is doable is that all the exogenous processes are known at time 0.

Basically, there are two balances, a and b, that get depleted over time. Say 1 starts at 800 and one starts at 200. Each period, there are two exogenous processes x and y. X depleted the balances a and b in proportion to their total balance. Y depletes y directly until it is 0 and then it depletes x.

For example, period 1: x = 10, and y = 5. End of period a = 800 - .8 * 10 = 792, y = 200 - .2 * 10 - 5 =193. The next period, the multipliers for x have changed. They will now be 792/985 for a and 193/995 for y.

Here is a very simple example using a loop:

data <- data.frame(start_a = NA, start_b = NA, proportion = NA, x = runif(10, 1, 50), y = runif(10, 1, 50), end_a = NA, end_b = NA)

for (i in 1:(nrow(data))){
    data$start_a[i] <- ifelse(i==1, 800, data$end_a[i-1])

    data$start_b[i] <- ifelse(i==1, 200, data$end_b[i-1])

    data$proportion[i] <- data$start_a[i]/(data$start_a[i]   data$start_b[i])

    data$end_a[i] <- data$start_a[i] - data$proportion[i]*data$x[i] -

    ifelse((data$start_b[i] - (1-data$proportion[i])*data$x[i])<= data$y[i], data$y[i] - data$start_b[i] - (1-data$proportion[i])*data$x[i], 0)

    data$end_b[i] <- data$start_b[i] - (1-data$proportion[i])*data$x[i] - min(data$y[i], data$start_b[i] - (1-data$proportion[i])*data$x[i])
}

Again, the full time history of x and y are known at the beginning, so I have a gut feeling the loop isn't needed.

CodePudding user response：

adding onto @jblood94 answer to use random vectors for x and y

n1 <- length(cumsum(y)[cumsum(y) < b])
aOut <- c(a b, a   b -cumsum(x y), 0)
aOut <- aOut[aOut > 0]
a1 <- c(a, a*cumprod(1 - x[1:n1]/aOut[1:n1]))
b1 <- aOut[1:length(a1)] - a1
idx <- match(TRUE, b1 < 0) - 1L
bOut <- c(b1[1:idx], rep(0, length(aOut) - idx))
aOut[1:idx] <- a1[1:idx]
output = data.frame(a = aOut, b = bOut)

CodePudding user response：

Here's a function that uses vectorization to get a and b according to your explanation:

fSeries <- function(a, b, x, y) {
  n1 <- b%/%y
  aOut <- seq(a   b, 0, -x - y)
  a1 <- c(a, a*cumprod(1 - x/aOut[1:n1]))
  b1 <- aOut[1:length(a1)] - a1
  idx <- match(TRUE, b1 < 0) - 1L
  bOut <- c(b1[1:idx], rep(0, length(aOut) - idx))
  aOut[1:idx] <- a1[1:idx]
  return(data.frame(a = aOut, b = bOut))
}

df <- fSeries(800, 200, 10, 5)
list(head = head(df), tail = tail(df))
#> $head
#>          a        b
#> 1 800.0000 200.0000
#> 2 792.0000 193.0000
#> 3 783.9594 186.0406
#> 4 775.8773 179.1227
#> 5 767.7530 172.2470
#> 6 759.5854 165.4146
#> 
#> $tail
#>     a b
#> 62 85 0
#> 63 70 0
#> 64 55 0
#> 65 40 0
#> 66 25 0
#> 67 10 0