#########################################
#  第14章 非ガウス型モデル(TSSSパッケージ)
#########################################
#
## trend model
#
#　テストデータの生成（n=101,201でジャンプ）
#
  x <- rep(0,400)
  x[101:200] <- 1
  x[201:300] <- -1
  y <- x + rnorm(400, mean=0, sd=1.0)
#
# system noise density : Gaussian (normal)
# 正規分布モデル：　カルマンフィルタと同じ結果
#
s1 <- ngsmth(y, noisev = 1, tau2 = 1.4e-02, noisew = 1, sigma2 = 1.048)

plot(s1, "smt", theta = 25, phi = 30, expand = 0.25, col="white")


#　システムノイズ：非ガウス分布
#　ピアソン分布(b=0.6)の場合
#
# system noise density : Pearson family
s2 <- ngsmth(y, noisev = 2, tau2 = 2.11e-10, bv = 0.6, noisew = 1,
             sigma2 = 1.042)

plot(s2, "smt", theta = 25, phi = 30, expand = 0.25, col="white")


#　コーシー分布の場合（b=1)
# system noise density : Pearson family
  ngsmth(y, 2, 3.53e-5,bv=1.0, 2, 1.045)
s3 <- ngsmth(y, noisev = 2, tau2 = 3.53e-5, bv = 1.0, noisew = 1,
             sigma2 = 1.045)

plot(s3, "smt", theta = 25, phi = 30, expand = 0.25, col="white")


##########################
#　時変分散モデル
#
##  an earthquake wave data
#
data(MYE1F)
n <- length(MYE1F)
yy <- rep(0, n)
for (i in 2:n) yy[i] <- MYE1F[i] - 0.5 * MYE1F[i-1]
m <- seq(1, n, by = 2)
y <- yy[m]
par(mar=c(2,2,1,1)+0.1)
z <- tvvar(y, trend.order = 2, tau2.ini = 4.909e-02, delta = 1.0e-06)

#
#
# system noise density : Gaussian (normal)
#　ガウス近似（カルマンフィルタと同じ結果）
#
s3 <- ngsmth(z$sm, noisev = 1, tau2 = z$tau2, noisew = 2, sigma2 = pi*pi/6,
             k = 190)

plot(s3, "smt", theta = 25, phi = 30, expand = 0.25, col = 8)