#########################################
#  第13章 時変分散と時変係数ARモデル
#########################################
#
# 地震波データの入力
#
data(MYE1F) # an earthquake wave data 

par(mar=c(2,2,1,1)+0.1)
# 
# 二重指数分布
#
ymean <- mean(MYE1F)
x <- MYE1F - ymean
y <- x^2
plot( y,type="l" )

z <- log10( y )
plot( z,type="l" )

#
# 二重指数分布(2個の二乗和)
#
N <- length(MYE1F)
N2 <- N/2
t <- 1:N2
y2 <- t
for (i in t) y2[i] <- (MYE1F[2*i-1]^2+MYE1F[2*i]^2)/2
plot(y2,type="l",lwd=2)

w2 <- log10(y2)
plot(w2,type="l",lwd=2)

###############################################
#　時変分散モデル (TSSSパッケージ)
#
# 時変分散モデルの推定とデータの等分散化
# トレンド次数=2，システムノイズ分散の初期値と探索幅
#
data(MYE1F) 

tvvar(MYE1F, 2, 6.6e-06, 1.0e-06)

# トレンド次数=1
tvvar(MYE1F, 1, 1.0e-01, 1.0e-02)



##############################################
#  時変係数ARモデル
#
#　地震データの入力
#
data(MYE1F) # an earthquake wave data 
#
#　時変係数ARモデルの推定
#　AR次数=4，AR係数トレンド次数=2，係数のジャンプなし
#

z <- tvar(MYE1F, trend.order = 2, ar.order = 8, span = 40,
           tau2.ini = 6.6e-06, delta = 1.0e-06)
z

#　時変スペクトル
par(mar=c(0,0,0,0)+0.1)  #作図枠の設定
spec <- tvspc(z$arcoef, z$sigma2)
plot(spec, theta = 30, phi = 40, expand = 0.5)

#  新しい3次元プロット（未公開）
#　plot.spec(spec)

#
#　n=630と1026の2か所で構造変化があったと仮定
#
z <- tvar(MYE1F, trend.order = 2, ar.order = 8, span = 40,
          outlier = c(630, 1026), tau2.ini = 6.6e-06, delta = 1.0e-06)
#
#　時変係数ARモデルから時変スペクトルを計算
#　時変スペクトルを3次元表示
#
z

spec <- tvspc(z$arcoef, z$sigma2)
plot(spec, theta = 30, phi = 40, expand = 0.5)



##########################################
#　時変スペクトルの3次元表示
##########################################
# seismic data
data(MYE1F)
z <- tvar(MYE1F, trend.order = 2, ar.order = 8, span = 20,
          outlier = c(630, 1026), tau2.ini = 6.6e-06, delta = 1.0e-06)
spec <-
tvspc(z$arcoef, z$sigma2)
#########################################
# 以下は試作品
 
# 最初のスペクトル(n=0)
nf <- 201
dt <- 2
dy <- 0.2
nf1 <- nf-dt
t <- 1:nf
tt <- 1:nf1
plot(t,spec$z[,1],type="l", xlim=c(1,501),ylim=c(-2,18)) y <- spec$z[,1] z
<- y
 
# 鳥瞰図(n=1,80)
nrep <- 1:80
for (i in nrep){
par(new=T)
t <- t+dt
  for (j in 1:nf) z[j] <- spec$z[j,i]+dy*(i-1)
  for (j in tt){
#  z[j] <- max(spec$z[j,i]+dy*(i-1),y[j+dt])
  z[j] <- max(z[j],y[j+dt])
  }
y <- z
plot(t,y,type="l", xlim=c(1,501),ylim=c(-2,18),xaxt='n',yaxt="n")