R のインストール
  http://cran.ism.ac.jp （統計数理研究所のミラーサイト）

本講義関連時系列解析用のRパッケージ

TSSS（Time Series analysis with State Space model）パッケージ
  http://jasp.ism.ac.jp/ism/TSSS/ 

TIMSACパッケージ
  http://jasp.ism.ac.jp/ism/timsac/ 

TTR パッケージ
options(repos="http://cran.ism.ac.jp")
install.packages("TTR")
library(TTR)

----------------------------------------------------------
####################
# データの読み込み #
####################

# 船舶のデータ
hakusan <- as.ts(read.csv("hakusan_new.csv"))
hakusan1 <- hakusan[,1]   # 方向角速度データ

# 太陽黒点数データ
sunspot <- as.ts(read.csv("sunspot_new.csv"))

# 東京の日最高気温データ
maxtemp <- as.ts(read.csv("maxtemp.csv"))

# アメリカの食品産業に従事する労働者人口
blsfood <- as.ts(read.csv("blsfood_new.csv"))

# 工業製品の卸売り高
whard <- as.ts(read.csv("whard_new.csv"))

# 地震データ（東西方向）
mye1f <- as.ts(read.csv("mye1f_new.csv"))

# 日経225データ
nikkei225 <- as.ts(read.csv("nikkei225_new2.csv",header=TRUE))

# 東京降雨データ（2年間の降雨日数）
rainfall <- as.ts(read.csv("rainfall_new.csv",header=TRUE))
rainfall2 <- rainfall[,4]/2  # 2年間の平均

# 榛原地下水位データ
haibara <- as.ts(read.csv("haibara_new.csv"))



#########################
### データのプロット　###
#########################
par(mar=c(2,2,1,1)+0.1) # 余白を小さくする
plot(hakusan1)
plot(sunspot)
plot(maxtemp,ylim=c(0,35))
plot(blsfood)
plot(whard)
plot(mye1f)
plot(nikkei225)
plot(lynx)
plot(haibara[,c(2,4)])   #榛原データの地下水位，気圧
plot(hakusan[,c(3,4,7)]) #船舶データの横揺れ，縦揺れ，舵角
plot(rainfall2,type="h") #2年間の降水割合


##############################
### WHARDデータの対数変換　###
##############################
plot(whard)
plot(log(whard))  # log(whard)のプロット


################################
### 日経225データの対数差分　###
################################
plot(diff(log(nikkei225)))
#　または
y <- log(nikkei225)
z <- diff(y)
plot(z)


##############################
### WHARDデータの季節階差　###
##############################
# 季節階差
y<-log(whard)
z<-rep(NA,n)
n <- length(whard)
period <- 12
z[1:n] <- NA
for(i in period+1:n){
z[i] <- y[i]-y[i-period]
}
plot(as.ts(z))

# 前期比
plot(whard/lag(whard))

# 前年同期比
plot(whard/lag(whard,k=-12))


################################
### 移動平均と移動メディアン ###
################################

# 移動平均フィルタ（TTRパッケージの関数 SMA)
plot(maxtemp,ylim=c(0,40))
x <- SMA(maxtemp,17)
lines(x,col=2,lwd=1)

# k=5,17,29の移動平均
plot(SMA(maxtemp,5),ylim=c(0,40))
plot(SMA(maxtemp,17),ylim=c(0,40))
plot(SMA(maxtemp,29),ylim=c(0,40))

# 移動平均フィルタ（自作版）

plot(maxtemp,ylim=c(0,40))
y <- maxtemp
ndata <- length(maxtemp)
y[1:ndata] <- NA
kfilter <- 17
n0 <- kfilter+1
n1 <- ndata-kfilter
for(i in n0:n1){
i0 <- i-kfilter
i1 <- i+kfilter
y[i] <- mean(maxtemp[i0:i1])
}
lines(y,col=2,ylim=c(0,40))

#
# 移動メディアンフィルタ
#
plot(maxtemp,ylim=c(0,40))
y <- maxtemp
ndata <- length(maxtemp)
y[1:ndata] <- NA
kfilter <- 17
n0 <- kfilter+1
n1 <- ndata-kfilter
for(i in n0:n1){
i0 <- i-kfilter
i1 <- i+kfilter
y[i] <- median(maxtemp[i0:i1])
}
lines(y,col=3,ylim=c(0,40),lwd=2)

#
# 移動平均・移動メディアン比較用テストデータ（ノイズなし・異常値あり）の生成
#
z<-rep(0,100)
z[51:100]<-1
z[24]<-1
z[75]<-2
z<-as.ts(z)
ploｔ（ｚ）
＃　移動平均
y <- z
ndata <- length(z)
y[1:ndata] <- NA
kfilter <- 3
n0 <- kfilter+1
n1 <- ndata-kfilter
for(i in n0:n1){
　i0 <- i-kfilter
　i1 <- i+kfilter
　y[i] <- mean(z[i0:i1])
}
lines(y,col=2,ylim=c(0,2),lwd=2)

# 移動メディアン
y <- z
ndata <- length(z)
y[1:ndata] <- NA
kfilter <- 3
n0 <- kfilter+1
n1 <- ndata-kfilter
for(i in n0:n1){
i0 <- i-kfilter
i1 <- i+kfilter
y[i] <- median(z[i0:i1])
}
lines(y,col=3,ylim=c(0,2),lwd=2)


#
# 移動平均・移動メディアン比較用テストデータの生成
#
x <- rep(0,400)
 x[101:200] <- 1
 x[201:300] <- -1
 y <- x + rnorm(400, mean=0, sd=0.5)
ng_test <-as.ts(y)
plot(ng_test)


#########################################
# 地下水位データの表示（異常値・欠測値）#
#########################################
# 榛原 地下水位データ
haibara <- as.ts(read.csv("haibara_new.csv"))
haibara_water <- haibara[,2]
plot(haibara_water,type="h")

# 一部分の詳細表示
plot(haibara_water[211:280],type="h",ylim=c(6.35,6.38))