####################################
#  第２章 共分散関数
####################################

###　データの読み込み
#
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"))

# 方向角データと横揺れデータ
plot(hakusan[,c(1,3)],main="")

###################################
# ヒストグラム
hist(hakusan[,1],main="方向角速度",col="#539952")
hist(hakusan[,3],main="横揺れ",breaks=seq(-6,12,1.2), col="#edae00")
hist(hakusan[,4],main="縦揺れ",breaks=seq(-20,20,3))
hist(hakusan[,7],main="舵角",breaks=seq(-15,5,1.5), col="#edae00")


###################################
#自己相関関数
unicor(hakusan1)   # TSSSパッケージの場合
unicor(subspot)
unicor(maxtemp)
unicor(blsfood)
unicor(whard)
unicor(mye1f)

acf(hakusan1)      #
acf(sunspot)
acf(maxtemp,lag.max=50)

#多変量時系列
#Package(TSSS)
crscor(hakusan[,c(3,4,7)])    #　船舶データ（横揺れ，縦揺れ，舵角）
crscor(hakusan[,c(1,3,4,7)])  #　方向角を加えた4変量

haibara <- read.csv("haibara_new.csv")  # 地下水位データ
haibara2 <- as.ts(haibara[,c(2,4)])     # 地下水位と気圧
crscor(haibara2)

# 多変量時系列の散布図
pairs(hakusan[,c(3,4,7)],col="blue")

hist(hakusan[,3],breaks=seq(-6,12,1.2),ylim=c(0,200),col="#539952")
hist(hakusan[,4],breaks=seq(-20,20,2.4),ylim=c(0,200),col="#539952")
hist(hakusan[,7],breaks=seq(-16,8,1.6),ylim=c(0,200),col="#539952")


##############################################
# 多変量正規乱数　　MASSパッケージ  mvrnorm
##############################################
# 独立ノイズ　r=0
Sigma <- matrix(c(1,0,0,1),2,2)
Sigma
x <- mvrnorm(n = 1000, rep(0, 2), Sigma)
plot(x, pch=19)

#　正の相関 r=0.6
Sigma <- matrix(c(1,-0.6,-0.6,1),2,2)
Sigma
x <- mvrnorm(n = 1000, rep(0, 2), Sigma)
plot(x, pch=19)

#####################
# AR process の生成
Sigma <- matrix(c(1,0,0,0,1,0,0,0,1),3,3)
Sigma

n <- 100
x <- mvrnorm(n = 100, rep(0, 3), Sigma)
y <- x
#x <- as.ts(x)
#plot(x)

for(j in 1:3}{
for(i in 3:n){
y[i,j] <- x[i,j] + 1.8*y[i-1,j] - 0.95*y[i-2,j]
}
}
y <- as.ts(y)
#####################


####################################
#  第３章 スペクトルとピリオドグラム
####################################

# periodogram
# hakusan data   船舶データ（方向角速度）
x <- period(hakusan1,window=0)
plot(x$smoothed,ylim=c(-3,2),type="n")
par(new=T)
y <- x$smooth+3
plot(y,ylim=c(0,5),type="h")

# sunspot data　　太陽黒点データ
x <- period(sunspot,window=0)
plot(x$smoothed,ylim=c(0,5),type="h")

# maxtemp data
x <- period(maxtemp,window=0)
plot(x$smoothed,ylim=c(0,4),type="h")

#blsfood
x <- period(blsfood,window=0)
plot(x$smoothed,ylim=c(0,5),type="h")

# whard
x <- period(whard,window=0)
plot(x$smoothed,ylim=c(0,7),type="h")

# mye1f
period(mye1f)

#############################
# FFT によるピリオドグラム
fftper(mye1f)
period(mye1f,plot=FALSE)

# 正規乱数の生成
r <- as.ts(rnorm(200))
plot(r)

# AR過程(2次）の生成
t <- 3:300
r <- as.ts(rnorm(300))
a1 <- 0.9*sqrt(3)
a2 <- -0.81
for (i in t)  x[i] = a1*x[i-1] + a2*x[i-2] + r[i]
plot(as.ts(x[101:300]))

# 正規白色雑音の生成(n=3200)
r <- as.ts(rnorm(3200))
plot(r)
acf(r)
plot(r,main="acf")

# ピリオドグラム（データ数3200）
r <- rnorm(3200)
r <- as.ts(r)
spec.pgram(r,log="yes")

# ピリオドグラム（データ数200）
spec.pgram(r[1:200],log="yes")

f <- 0:113
f <- f/226
x <- as.ts( period(r[1:200]) )
plot(f,x,type="l")


# ピリオドグラム（データ数800）
#spectrum
spectrum(r,ylim=c(0.01,100))
spectrum(r[1:800],ylim=c(0.01,100))


##  Raw Spectrum
# autocovariance function
lag <- 50
lag <- length(r)
lag1 <- lag-1
x <- as.ts(acf(r,type="covariance",lag.max=lag))
c <- x$acf

# raw spectrum
p <- rep(1:lag)
f <- 0:lag
g <- 1:lag1
for (i in f)  {
  fj <- i/(2*lag)
  s <- c[1]
  for (k in g) {
     s <- s + 2*c[k+1]*cos(2*pi*k*fj)
     }
  p[i+1] <- s
}
p <- as.ts(p)
lp <- log10(p)
plot(lp,ylim=c(-4,1))