重新运行公式以生成许多列数据（矩阵使用函数rep（）＆Matrix（））

我最初使用的R代码低于其中使用公式中的常量来运行模型，并使用一个常量数字wb $ yfit和X创建数据模拟，其中每个数据点都包含一个随机正态变量 . 目标是不仅使用创建一列数据，而且使用Z行的列数，其中行与wb $ yfit重合，但由于每个[，＃]处的X而变化 .

N <- length(wb$time) # No. of historical observations
theta <- theta.hat
sigma <- sigma.hat
stdev <- sqrt ( sigma ^2/(2* theta )*(1- exp (-2* theta )))
X <- vector(mode="numeric", length =N)
X[1] <- wb$Ln_price_[1] - wb$yfit[1]

for (i in 2:N){
    X[i] <- exp (- theta ) * X[i -1] + stdev * rnorm (1 ,0 ,1)
}

NBP.sim <- wb$yfit + X
plot (wb$time, NBP.sim, ylim =c (3 ,5) ,type ="l",
      col="red",xlab =" Time ",ylab ="EUR / MWh ",las =1)
lines (NBP.sim , col=" blue ", type ="l")

编码它的想法是这样的，但我不能让它工作：

N <- 10
days <- length(wb$time)
X <- rep(wb$Ln_price_[1] - wb$yfit[1],N)
TTF.sim.Q <- matrix(0, nrow =N, ncol=days)

for(i in 2: days ){
    ## Simulating the OU part of the log price
    X[i] <- exp (- theta ) * X[i-1] + stdev * rnorm (N ,0 ,1)
    ## Adding the deterministic part to the stochastic 
    ## part to get the final simulated spot price
    TTF.sim.Q[,i] <- wb$yfit[i] + X 
}

## Plotting prices in EUR /MWh
matlines(wb$time, t(TTF.sim.Q),type ="l",
         lty= rep (1, ncol(t(TTF.sim.Q))), ylim =c (3 ,5))

非常感谢!!!

我最初使用的R代码使用公式来运行模型，并使用数字wb $ yfit和X创建数据模拟，其中涉及每个数据点的随机正态变量 . 这导致一列731模拟价格，使用wb $ yfit并添加X的随机流量 . 这个X是标准偏差和价格的均值回归成分，因此随机变量的原因是价格上涨或下跌随意的 .

*代码的最后一部分在一个图上绘制了这个模拟，以显示价格随时间的变化 .

目标是不仅使用创建一列数据，而且使用行数（时间）与wb $ yfit（不同时间的不同值）重合的Z列数量，但由于X包含随机变量而因X而变化 .

那么我怎样才能进行价格的10次模拟，而不仅仅是我创建的模拟，这在我看来意味着我想创建10个模拟，意味着使用上面代码中使用的相同函数的10列数据 . 每列代表731价格随时间的模拟（每个点为wb $ yfit [n] X [n]，下一个点将是wb $ yfit [n 1] X [n 1]） . 因此，由于计算中包含的随机变量，列中每个时间点的价格应该不同 . X.最后我想在图中显示这10个模拟 .

知道我怎么能创造这个吗？

我假设它仍然使用rep（）和Matrix（） . （我想最大的问题是为每个X生成一个独立的随机＃，取决于之前的X，以及将其添加到相关WB $ yfit的时间）再次感谢..

N <- length(wb$time) # No. of historical observations
theta <- theta.hat
sigma <- sigma.hat
stdev <- sqrt ( sigma ^2/(2* theta )*(1- exp (-2* theta )))
X <- vector(mode="numeric", length =N)
X[1] <- wb$Ln_price_[1] - wb$yfit[1]

for (i in 2:N){
    X[i] <- exp (- theta ) * X[i -1] + stdev * rnorm (1 ,0 ,1)
}

NBP.sim <- wb$yfit + X
plot (wb$time, NBP.sim, ylim =c (3 ,5) ,type ="l",
      col="red",xlab =" Time ",ylab ="EUR / MWh ",las =1)
lines (NBP.sim , col=" blue ", type ="l")

编码它的想法是这样的，但我不能让它工作：

N <- 10
days <- length(wb$time)
X <- rep(wb$Ln_price_[1] - wb$yfit[1],N)
TTF.sim.Q <- matrix(0, nrow =N, ncol=days)

for(i in 2: days ){
    ## Simulating the OU part of the log price
    X[i] <- exp (- theta ) * X[i-1] + stdev * rnorm (N ,0 ,1)
    ## Adding the deterministic part to the stochastic 
    ## part to get the final simulated spot price
    TTF.sim.Q[,i] <- wb$yfit[i] + X 
}

## Plotting prices in EUR /MWh
matlines(wb$time, t(TTF.sim.Q),type ="l",
         lty= rep (1, ncol(t(TTF.sim.Q))), ylim =c (3 ,5))

非常感谢!!!

我到现在为止尝试的是：

N < - 10

天< - 长度（wb $时间）X < - rep（wb $ Ln_price_ [1] - wb $ yfit [1]，N）TTF.sim.Q < - 矩阵（0，nrow = N，ncol =天） for（i in 2：days）{X [i] < - exp（-theta）* X [i-1] stdev * rnorm（N，0,1）＃模拟原木价格的OU部分

TTF.sim.Q [，i] < - wb $ yfit [i] X＃将确定性部分添加到随机部分以获得最终模拟现货价格}

TTF.sim.Q [，i]中的错误< - wb $ yfit [i] X：要替换的项目数不是替换长度的倍数另外：警告消息：1：在X [i] < - exp（ -theta）* X [i - 1] stdev * rnorm（N，0,1）：要替换的项目数不是替换长度的倍数2：在X [i] < - exp（-theta）* X [ i - 1] stdev * rnorm（N，0,1）：要替换的项目数不是替换长度的倍数3：在X [i] < - exp（-theta）* X [i - 1] stdev * rnorm（N，0,1）：要替换的项目数不是替换长度的倍数4：在X [i] < - exp（-theta）* X [i - 1] stdev * rnorm（N，0， 1）：要替换的项目数不是替换长度的倍数5：在X [i] < - exp（-theta）* X [i - 1] stdev * rnorm（N，0,1）：项目数替换不是替换长度的倍数6：在X [i] < - exp（-theta）* X [i - 1] stdev * rnorm（N，0,1）：要替换的项目数不是多个替换长度7：在X [i] < - exp（-theta）* X [i - 1] stdev * rnorm（N，0,1）：数字o f要替换的项目不是替换长度的倍数8：在X [i] < - exp（-theta）* X [i - 1] stdev * rnorm（N，0,1）：要替换的项目数不是替换长度的倍数9：在X [i] < - exp（-theta）* X [i - 1] stdev * rnorm（N，0,1）：要替换的项目数不是替换长度的倍数10 ：在X [i] < - exp（-theta）* X [i - 1] stdev * rnorm（N，0,1）：要替换的项目数不是替换长度的倍数

检查结果 - >

data.frame（X）X 1 -0.07735493 2 -0.08452771 3 -0.09945222 4 -0.04269141 5 -0.03170452 6 -0.04693956 7 0.01079676 8 0.03830161 9 0.06544999 10 0.08219423 11 0.09241736

当我尝试代码时：

N < - 10

天< - 长度（wb $时间）X < - rep（wb $ Ln_price_ [1] - wb $ yfit [1]，N）TTF.sim.Q < - 矩阵（0，nrow = N，ncol = days）erri < - rnorm（length（days））for（i in 2：days）{##模拟原木价格的OU部分X [i] < - exp（ - theta）* X [i-1 ] STDEV * ERRI [I] ##部得到最终的模拟现货价格TTF.sim.Q [I] < - WB $ yfit [I] X [I]} data.frame（X）X 1 -0.07735493 2 NA 3 NA

和NA一样......

TTF.sim.Q [i] [1] NA

而TTF.sim.Q什么也没给我