Ошибка в петле в R, которую я не могу найти

Question

У меня есть следующая модель, которую я хочу рассчитать: где

Мои данные следующие:

D0= 100
variance = 3
T = 20
uf = 0.3
ue = 0.7
CRRA = 0.1
theta = 0.7
Q = 1
X1=1

# DIVIDENDS:
epsilon <- c(0,0, 6, rep(0,18))

DD_t <- D0
for (t in 2:(T+1)) {
    DD_t[t] <- epsilon[t]+DD_t[t-1]
}

# PRICE in t0 t1 t2
PP_t<- c(rep(0,21))
PP_t[1] <- DD_t[1] - 0.1*(3^2)*(18 + (1/0.3))
PP_t[2] <- DD_t[2] + (ue/uf)*theta^(1-1) - CRRA*(variance^2)*Q*(T-1-1 + (1/uf))
PP_t[3] <- DD_t[3] + (ue/uf)*theta^(2-1) - CRRA*(variance^2)*Q*(T-2-1 + (1/uf))

До этого момента все правильно. Теперь я хочу вычислить от P_3 до P_20, и это следующий код, к сожалению, дает неправильное решение.

# PRICE in t3 t4 to t20

for (t in 3:(T-1)) {
    t1 <- (t + 1)

    XX_t <- 0
    for (k in 1:(t-1)) {
    XX_t <- XX_t + (theta^(k-1))*(PP_t[t1-k]-PP_t[t1-k-1]) 
    }

    XX_t <- (1-theta)*XX_t + theta^(t-1)*1
    PP_t[t1] <- DD_t[t1] + XX_t*(ue/uf) - CRRA*(variance^2)*Q*(T - t - 1 + (1/uf))
}

Что это дает мне:

PP_t

 [1]  80.80000  83.13333  89.33333  95.22667  98.26400  98.66093  97.36051  95.53206  94.15252  93.77109  94.47277  95.99211  97.89801  99.77274 101.33323 102.47560 103.25258
[18] 103.81069 104.31816 104.90794   0.00000

Однако правильное решение выглядит так:

 [1]  80.80000    83.13333  89.33333  94.08333  96.66333  97.22033  96.46413  95.28555  94.43488  94.33895  95.07011  96.43089  98.09417  99.74021 101.15038 102.24238 103.05407
[18] 103.69603 104.29659 104.95801 105.73360

К сожалению, я не могу найти ошибку в своем коде. Оно должно быть незначительным, поскольку значения не так уж и далеки.