eta <- rorm(T, 0, sqrt(tau2)) # AR(1)波动率模型的误差
    z <- rnrm(T, 0, 1) # 倍增项回报模型
    alha\[1\] <- cost # 在开始阶段没有自相关的观察值

# 仿真时间序列
smdf <- s_sm(theta, T)
y <- smdf$y
lpa <- smdf$apha

隐马尔可夫模型：定义

上面显示的模型属于更一般的隐马尔可夫模型类。设 h(αt|αt-1;θ) 为跃迁密度，g(yt|αt;θ) 为测量密度。那么在这种情况下，跃迁密度和测量密度都是高斯的，其中  和 .

序列蒙特卡罗

对于估计，我们使用序列蒙特卡罗，通过生成 P 随机抽取，称为“粒子”，以近似预测和过滤密度。虽然有很多变体，但我们只讨论（连续）序列重要性重采样（SIR）。

具有连续序列重要性重采样的过滤步骤：算法

连续序列重要性重采样(CSIR) 是 SIR 的一种变体，它提供了过滤粒子的连续版本。

使用 CSIR 的过滤步骤的算法如下：

• 输入：
• 具有条目 u(j) 的排序均匀随机采样向量（拒绝采样）；
• 对于定义为 W(i)t 的每个粒子 α(i)t 在 yt 处评估的正态 PDF；
• 从预测密度 α(i)t 中排序。

代码

下面我们生成粒子集，并使用 SIR 近似过滤和预测密度。在第一个图中，我们显示了预测密度平均值及其 95 和 5 分位数。在同一个图中，我们还绘制了波动率的真实值。在第二个图中，我们绘制了过滤密度的热图。黑线是真正的波动率。

# --> （原始）序列重要性取样算法：过滤步骤

# 输入 1: appr - 预测密度
# 输入 2: aha_t - 在 y\[t\]评估的正态 pdf
# 输入 3: u - 排序均匀的随机采样向量(拒绝采样)
# 输出：alphp - 粒子过滤


    # 排序和加权的速度减慢
    alhawt <- alph\_wt/sum(alpha\_wt)
    alpa\_rt <- cbind(seq(1,P,1),alpha\_pr)
    alhapr\_id <- lpha\_sort\[order(alha_r\[,2\]),\]。
    alhapr <- alpha_ridx\[,2\]
    alph\_ix <- alha\_p_idx\[,1\] 
    alha\_wt <- alp\_w\[alpha_idx\]
    alhacwt <- c(0, cumsum(alpha_wt))

    j <- 1
    for (i in 1:P)
        while((aphawt\[i\] < u\[j\]) && (u\[j\] <= alpawt\[i+1\])){
            lp\_up\[j\] <- alpa\_r\[i\] 。


}

# ----------------------------------------------------------------------
# 设置粒子过滤
# ----------------------------------------------------------------------
P <- 200 # 设置粒子的数量
lph_up <- rnorm(P,0,0.1)
alpar <- rep(0,P)
aha_w <- rep(1,P)/P

alphup_mt <- matrix(rep(0,T*3),T)
ala_pmat <- matrix(rep(0, T*3),T)
ah_prare <- matrix(rep(0, T*20),T)


# 从一个近似值中生成一个由P个随机抽样组成的粒子集
# 每个时间序列点的预测和过滤分布的近似值
for (t in 1:T){
    # 预测步骤
    appr <- nst + phi * alpp + rnorm(P,0,srt(tau2))
    # 更新/过滤步骤(态密度)
    ahat <- dnorm(y\[t\]*rep1,P), mean=0 , sd = exp(phar/2)
    alpap <- sir(alhapr=aph\_r,alhawt=alpa\_t, u=sort(runif(P,0,1))
# 绘制预测密度图

plot(sqrt(252) * exp(alpha/2), type='l')

## 筛选密度热图
    
het <- matrix(rep(1,T*20), T, 20)

plot(NULL, xlim = c(1, T), ylim = c(0, 160), main="过滤密度热图",

###连续序列重要性重取样：过滤步骤

# 输入 1: alppr - 预测密度
# 输入 2: alhawt - 在 y\[t\]处评估的正态 pdf
# 输入 3: u - 排序均匀的随机采样向量(拒绝采样)
# 输出：ala_up - 粒子过滤（连续版本）。

# R版本(性能较慢)
cir <- function(aph_r, phwt, u) {
    P <- length(aphpr)

    al_p <- rep(0,P)
    
    # 排序和加权的速度减慢
    alpha\_wt <- alpha\_wt/sum(alpha_wt)
  

    j <- 1
    for (i in 1:P){
        while((a_ct\[i\] < u\[j\]) & (u\[j\] <= alhwt\[i+1\])){
            alh\_u\[j\] <- aph\_pr\[i\] + ((apapr\[i+1\]-alar\[i\])/(ala\_ct\[i+1\]-alpha\_cwt\[i\]) * (u\[j\]-ala_wt\[i\])
         
}


csir.c <- function(alppr, aht, u) {
    P <- length(alpap)
    ala_u <- rep(0,P)
     .C("cir", alpup=as.dole(aphup), 
                          alha\_pr=as.double(aha\_r), 
                          alh_wt=as.doublephawt), 
                          u=as.double(u),

我们现在提供用于最大化对数似然和估计参数 θ 的代码。为了计算标准误差，我们使用在 MLE 评估的对数似然的 Hessian 矩阵的逆矩阵的对角线。

我们现在可以转到参数 θ 的估计。使用 C 中的函数进行估计。

vas <- sfit(y, c(0.5,0.5,0.5), P, 1)

## 显示结果
matrix <- cbind(heta_mle
,eta_se)
矩阵