#### 第6章 ####
#### 6.1節 ####
# パネルデータの時系列分析
# パネルデータ作成
T=50 # 期間
n=3  # 個別主体数
rho=0.8
sig1=1.2
y=matrix(0,nrow=T,ncol=n)
X=matrix(0,nrow=T,ncol=n)
X[1,]=rnorm(n,0,1)
for(time in 2:T){
for(i in 1:n){
X[time,i]=rho*X[time-1,i]+rnorm(1,0, 1/sqrt(sig1))}}
time<-1:T
plot(time,X[,1], type = "l", main = "Panel data", xlab = "time", ylab = "x(t)")
for(i in 2:n){lines(time,X[,i],col=i)}

sig2<-1
beta=matrix(rnorm(n,0,1/sqrt(sig2)),nrow=n,ncol=1)
for(i in 1:n){
y[,i]<-X[,i]*beta[i]+rnorm(50,0,1/sqrt(sig2))}
panel.lm <- lm(c(y)~c(X)-1)
#panel.lm <- lm(y~X-1)
summary(panel.lm)

# データのエクスポート
colnames(X)=c("x1", "x2","x3")
rownames(X)=c(time)
colnames(y)=c("y1","y2","y3")
rownames(y)=c(time)
write.table(X,"C:/Program Files/R/R-2.7.2/panelX.csv",sep=",")
write.table(y,"C:/Program Files/R/R-2.7.2/panelY.csv",sep=",")

#### 個別主体数を３に固定する場合####
x1<-x2<-x3<-real(length=n)
x1[1]<-0.7; x2[1]<--0.2; x3[1]<-0.3
for(time in 2:T){
x1[time]<- rho*x1[time-1] + rnorm(1, 0, 1/sqrt(sig1))
x2[time]<- rho*x2[time-1] + rnorm(1, 0, 1/sqrt(sig1))
x3[time]<- rho*x3[time-1] + rnorm(1, 0, 1/sqrt(sig1))}
X=cbind(x1,x2,x3)
time<-1:T
plot(time,x1, type = "l", main = "AR(1) process", xlab = "time", ylab = "x(t)")
lines(time,x2,col=2)
lines(time,x3,col=3)
beta1<-0.3; beta2<-1.2; beta3<--0.6; sig2<-1
y<-beta1*x1+beta2*x2+beta3*x3+rnorm(50,0,1/sqrt(sig2))
panel.lm <- lm(y~X-1)
summary(panel.lm)

# BUGSコード by Lancaster(2006) "panel.txt"を加筆修正
model{
for(i in 1:n){for(tt in 1:time){
y[tt,i]~dnorm(mu[tt,i],tau)
mu[tt,i]<-beta*X[tt,i]+alpha[i]
alpha[i]~dnorm(alphabar,phi)}}
for(i in 1:n){alpha[i]~dnorm(alphabar,phi)}
alphabar~dnorm(0,0.0001)
beta~dnorm(0,0.0001)
tau~dgamma(0.01,0.01)
phi~dgamma(0.01,0.01)}

# パネルデータ分析
library(R2WinBUGS)
X=as.matrix(read.table("C:/Program Files/R/R-2.7.2/panelX.csv",sep=",",header=T))
y=as.matrix(read.table("C:/Program Files/R/R-2.7.2/panelY.csv",sep=",",header=T))
time=nrow(X)
n=ncol(X)

data<-list("n","time","y","X")

inits<-function(){
list(
alpha=rnorm(n,0,0.0001),
alphabar=rnorm(1,0,0.0001),
phi=rgamma(1,0.1,0.1),
beta=rnorm(1,0,0.1),
tau=rgamma(1,0.1,0.1))}

parameters<-c("alpha","alphabar","beta","tau","phi")

panel.post<-bugs(data,inits,parameters,
model.file="C:/Program Files/R/R-2.7.2/panel.txt",debug=FALSE,
n.iter=10000,n.burnin=1000,n.chains=5,
bugs.directory = "C:/Program Files/WinBUGS14",
working.directory = "C:/Program Files/R/R-2.7.2")
panel.post
# 図6.1
plot(panel.post)

#### 6.2節 ####
# AR(1)
# データ作成
# Lancaster(2006)のコードを参照
T<-50
rho<-0.8
sig<-1.2
y<-real(length=T)
y[1]<-rnorm(1,0,0.1)
for(i in 2:T){y[i]=rho*y[i-1]+rnorm(1,0,1/sqrt(sig))}
time<-1:T
ar1=as.matrix(cbind(time,y))
write.table(ar1,"ar1.csv",sep=",")
rm(ar1)
ar1<-read.table("ar1.csv",sep=",",header=T)
time=ar1$time
y=ar1$y
plot(time,y, type = "l", main = "AR(1) process", xlab = "time", ylab = "y(t)")

# BUGSコード 
# ar1_1.txt Lancaster(2006)のコードを参照
model{
for(t in 2:T){
y[t]~dnorm(m[t],tau)
m[t]<-mu*(1-rho)+rho*y[t-1]}
rho~dunif(-.999,.999)
mu~dnorm(0,.001)
tau~dgamma(.001,.001)}

# ar1_2.txt
model{
for(t in 2:T){
y[t]~dnorm(m[t],tau)
m[t]<-alpha+rho*y[t-1]}
alpha~dnorm(0,.001)
rho~dnorm(0,.001)
tau~dgamma(.001,.001)
}

# Rコード
library(R2WinBUGS)
ar1<-read.table("C:/Program Files/R/R-2.7.2/ar1.csv",sep=",",header=T)
time=ar1$time
y=ar1$y
# 図6.2
plot(time,y, type = "l", main = "AR(1) process", xlab = "time", ylab = "y(t)")

T=nrow(ar1)
data<-list("T","y")

inits<-function(){
list(
alpha=rnorm(1,0,0.1),
rho=rnorm(1,0,0.1),
tau=rgamma(1,0.1,0.1))}

parameters<-c("alpha","rho","tau")

ar1.post<-bugs(data,inits,parameters,
model.file="C:/Program Files/R/R-2.7.2/ar1_2.txt",debug=FALSE,
n.iter=10000,n.burnin=1000,n.chains=3,
bugs.directory = "C:/Program Files/WinBUGS14",
working.directory = "C:/Program Files/R/R-2.7.2")
ar1.post
# 図6.3
plot(ar1.post)

# ar2.txt 予測こみ Congdon(2006)のコードを参照
# 動作未確認
model{
for (t in 1:T) {
y[t]~dnorm(mu[t],tau)}
mu[1]<-beta[1]+M[2]
mu[2]<-beta[1]+beta[2]*y[1]+M[1]
for(t in 3:T){
mu[t]<-beta[1]+beta[2]*y[t-1]+beta[3]*y[t-2]}
for (j in 1:4){
y[T+j]~dnorm(mu.new[T+j],tau)
mu.new[T+j]<-beta[1]+beta[2]*y[T+j-1]+beta[3]*y[T+j-2]}
for(j in 1:3){
beta[j]~dnorm(0.0,1.0E-2);}
for(j in 1:2){
M[j]~dnorm(0.0,1.0E-2);}
tau~dgamma(0.1,0.1)}

# Rコード
library(R2WinBUGS)
T=50
y=c(y,"NA","NA","NA","NA")
data<-list("T","y")

inits<-function(){
list(
M=c(0,0),
beta=c(4,0,0))}

parameters<-c("M","beta")

panel.post<-bugs(data,inits,parameters,
model.file="ar2.txt",debug=FALSE,
n.iter=10000,n.burnin=1000,n.chains=5,
bugs.directory = "C:/Program Files/WinBUGS14",
working.directory = NULL)
panel.post

time2=1:100
plot(1:100,y, type = "l", main = "AR(1) process", xlab = "time", ylab = "y(t)")
lines(time2,yp2,col=2)

#### 6.3節 ####
#ARMA(1,1) ARMA1_1.txt
model{
for(t in 1:T){y[t]~dnorm(m[t],tau.e)}
m[1]<-u[1]+mu
for(t in 2:T){m[t]<-rho*y[t-1]+u[t]+thita*u[t-1]}
for(t in 1:T){u[t]~dnorm(0,tau.u)}
mu~dnorm(0,0.001)
rho~dnorm(0,0.001)
thita~dnorm(0,0.001)
tau.e~dgamma(.001,0.001)
tau.u~dgamma(.001,0.001)
}

#Rコード
library(R2WinBUGS)
arma<-read.table("C:/Program Files/R/R-2.7.2/ar1.csv",sep=",",header=T)
time=arma$time
y=arma$y
T=nrow(arma)
data<-list("T","y")

inits<-function(){
list(
mu=rnorm(1,0,0.1),
rho=rnorm(1,0,0.1),
thita=rnorm(1,0,0.1),
tau.e=rgamma(1,0.1,0.1),
tau.u=rgamma(1,0.1,0.1))}

parameters<-c("mu","rho","thita","tau.e","tau.u")

arma.post<-bugs(data,inits,parameters,
model.file="C:/Program Files/R/R-2.7.2/arma1_1.txt",debug=FALSE,
n.iter=10000,n.burnin=1000,n.chains=3,
bugs.directory = "C:/Program Files/WinBUGS14",
working.directory = "C:/Program Files/R/R-2.7.2")
arma.post
# 図6.4
plot(arma.post)

#### 6.4節 ####
# VAR(4)モデル "var_4.txt" (arranged by congdon(2006) program8.5)
model{
for(t in 1:T) {
C[t] <- CONS[t]
Y[t] <- INC[t]}

for(t in p+1:T){
C[t] ~ dmnorm(C.mu[t], C.tau)
Y[t] ~ dmnorm(Y.mu[t], Y.tau)
C.mu[t] <- alp1+Clag1[t]+Ylag1[t]
Y.mu[t] <- alp2+Ylag2[t]+Clag2[t]
Clag1[t] <-  rho1[1]*C[t-1]+ rho1[2]*C[t-2]+ rho1[3]*C[t-3]+ rho1[4]*C[t-4]
Ylag1[t] <- beta1[1]*Y[t-1]+beta1[2]*Y[t-2]+beta1[3]*Y[t-3]+beta1[4]*Y[t-4]
Ylag2[t] <-  rho2[1]*Y[t-1]+ rho2[2]*Y[t-2]+ rho2[3]*Y[t-3]+ rho2[4]*Y[t-4]
Clag2[t] <- beta2[1]*C[t-1]+beta2[2]*C[t-2]+beta2[3]*C[t-3]+beta2[4]*C[t-4]}

alp1 ~ dnorm(0,0.001)
alp2 ~ dnorm(0,0.001)

for(m in 1:p){
rho1[m] ~ dnorm(0,0.01)
rho2[m] ~ dnorm(0,0.01)
beta1[m] ~ dnorm(0,0.01)
beta2[m] ~ dnorm(0,0.01)}

C.tau ~ dgamma(0.001,0.001)
Y.tau ~ dgamma(0.001,0.001)}

# Rコード
library(R2WinBUGS)
library(tseries)
jpyc<-read.table("C:/Program Files/R/R-2.7.2/jpyc.csv",sep=",",header=T)
jpyc.d=matrix(0,nrow=(nrow(jpyc)-5),ncol=2)
for(i in 1:(nrow(jpyc)-5)){
jpyc.d[i,1]<-log(jpyc[i+4,2])-log(jpyc[i+3,2])
jpyc.d[i,2]<-log(jpyc[i+4,3])-log(jpyc[i+3,3])
}
jpyc.ts<-ts(jpyc.d,start=c(1995),frequency=4)
colnames(jpyc.ts) <- c("GDP","CONS")
ts.plot(jpyc.ts,col=c(1,2))

# VAR(4)の最尤推定
jpyc.var<-ar(jpyc.ts,aic=T,order.max=4)
plot(jpyc.var$res)

# VAR(4)のベイズ推定
INC=c(jpyc.ts[,1])
CONS=c(jpyc.ts[,2])
T=length(INC)
p=4 # ラグ

data<-list("T","INC","CONS","p")
inits<-function(){
list(
alp1=rnorm(1,0,1.0),
alp2=rnorm(1,0,1.0),
rho1=rnorm(p,0,1.0),
rho2=rnorm(p,0,1.0),
beta1=rnorm(p,0,1.0),
beta2=rnorm(p,0,1.0),
C.tau=rgamma(1,0.1,0.1),
Y.tau=rgamma(1,0.1,0.1))}

parameters<-c("alp1","alp2","rho1","rho2","beta1","beta2","C.tau","Y.tau")
var4.post<-bugs(data,inits,parameters,
model.file="C:/Program Files/R/R-2.7.2/var_4.txt",debug=FALSE,
n.iter=10000,n.burnin=1000,n.chains=3,
bugs.directory = "C:/Program Files/WinBUGS14",
working.directory = "C:/Program Files/R/R-2.7.2")
var4.post
# 図6.7
plot(var4.post)

#### 6.5節 ####
# ARCH(1)/GARCH(1,1)
library(tseries)
library(fSeries)
# arch(1)およびgarch(1,1)分析用データの生成
n <- 600
a <- c(0.2,0.5) 
e <- rnorm(n)
y <- double(n)
y[1] <- rnorm(1, sd = sqrt(a[1]/(1.0-a[2])))
for(i in 2:n){
y[i] <- e[i]*sqrt(a[1]+a[2]*y[i-1]^2)}
y <- ts(y[101:n])

# ARCH(1)の最尤推定
y.arch <- garch(y, order = c(0,1)) 
summary(y.arch) 
# 図6.8
plot(y.arch)

# ARCH(1) arch_1.txt
# BUGSコード (Congdon (2006)を加筆修正) 
# テキストではOpenBugsを使っているが、WinBUGSを使ってもよい
# 以下では、WinBUGSを使った方法を示す
# OpenBUGSは、以下のサイトからダウンロードして利用できる
# http://mathstat.helsinki.fi/openbugs/

# 期間Tをあらかじめ指定しておく
model {for (t in 2:500) {
y[t] ~ dnorm(0,P[t])
h[t] <- gam+alph*pow(y[t-1],2)
P[t] <- 1/h[t]}
gam ~ dgamma(1,1)
alph ~ dunif(0,1)}

# Rコード
library(R2WinBUGS)
y=c(y)
data<-list("y")
inits<-function(){
list(
gam = rgamma(1,1,1),
alph = runif(1,0,1))}
parameters<-c("gam","alph")
arch1.post<-bugs(data,inits,parameters,
model.file="C:/Program Files/R/R-2.7.2/arch_1.txt",
n.iter=10000,n.burnin=1000,n.chains=2,
#program = "openbugs",
bugs.directory = "C:/Program Files/WinBUGS14/",
#bugs.directory = "C:/Program Files/OpenBUGS/",
working.directory = "C:/Program Files/R/R-2.7.2", digits=10)
arch1.post
# 図6.9
plot(arch1.post)

# GARCH(1,1)の最尤推定
y.garch <- garch(y, order = c(1,1)) 
summary(y.garch) 
# 図6.10
plot(y.garch)

# GARCH(1,1)（Congdon (2006)を参照）
# 誤差項が正規分布のモデル "garch1_1n.txt"
model {for (t in 2:500) {
y[t] ~ dnorm(0,P[t])
h[t] <- gam+alph*pow(y[t-1],2)+beta*h[t-1]
P[t] <- 1/h[t]}
h[1] <- h.1
h.1 ~ dgamma(1,1)
gam ~ dgamma(1,1)
alph ~ dnorm(0,1) I(0,)
beta ~ dnorm(0,1) I(0,)
sumcoef <- alph+beta
Stat <- step(1-sumcoef)}

# Rコード
library(R2WinBUGS)
set.seed(66)
y=c(y)
data<-list("y")
inits<-function(){
list(
gam = rgamma(1,.1,.1),
alph = runif(1,0,.01),
beta = runif(1,0,.01),
h.1 = rgamma(1,.1,.1))}
parameters<-c("gam","alph","beta","h.1")
garch1n.post<-bugs(data,inits,parameters,
model.file="C:/Program Files/R/R-2.7.2/garch1_1n.txt",
n.iter=10000,n.burnin=1000,n.chains=2,
program = "openbugs",
bugs.directory = "C:/Program Files/WinBUGS14/",
#bugs.directory = "C:/Program Files/OpenBUGS/",
working.directory = "C:/Program Files/R/R-2.7.2", 
digits=5)
garch1n.post
plot(garch1n.post)

# 誤差項がｔ分布のモデル "garch1_1t.txt"
model {for (t in 2:500) {
y[t] ~ dt(0,P[t],k)
h[t] <- gam+alph*pow(y[t-1],2)+beta*h[t-1]
P[t] <- 1/h[t]}
h[1] <- h.1
h.1 ~ dgamma(1,1)
gam ~ dgamma(1,1)
alph ~ dnorm(0,1) I(0,)
beta ~ dnorm(0,1) I(0,)
sumcoef <- alph+beta
Stat <- step(1-sumcoef)
k ~ dexp(0.1)}

# Rコード
library(R2WinBUGS)
set.seed(66)
y=c(y)
data<-list("y")
inits<-function(){
list(
gam = rgamma(1,0.1,0.1),
alph = runif(1,0,0.01),
beta = runif(1,0,0.01),
h.1 = rgamma(1,0.1,0.1),
k = rexp(1,0.1))}
parameters<-c("gam","alph","beta","h.1","k")
garch1t.post<-bugs(data,inits,parameters,
model.file="C:/Program Files/R/R-2.7.2/garch1_1t.txt",
n.iter=1000,n.burnin=100,n.chains=2,
#program = "openbugs",
bugs.directory = "C:/Program Files/WinBUGS14/",
#bugs.directory = "c:/Program Files/OpenBUGS/",
working.directory = "C:/Program Files/R/R-2.7.2", 
digits=5)
garch1t.post
# 図6.11
plot(garch1t.post)

# SVモデル(Lancaster (2006)を参照)
# データ作成
T=1000;mu1=0; rho=0.5; sig=1; h=rep(0,T)
h[1]=rnorm(1,mu1/(1-rho),1/sqrt(sig*(1-rho^2)))
for(j in 2:T){h[j]=rnorm(1,mu1+rho*h[j-1],1/sqrt(sig))}
y=exp(h/2)*rnorm(T)
time=1:T
svdat=cbind(time,y,h)
write.table(svdat,"C:/Program Files/R/R-2.7.2/svdat.csv",sep=",")
rm(svdat)
plot(time,h, type = "l", main = "Stochastic Volatility", xlab = "time", ylab = "SV")

# BUGSコード
model{
for(i in 1:T){
y[i]~dnorm(0,p[i])
p[i]<-exp(-h[i])}
h[1]~dnorm(muh,qh)
muh<-mu1/(1-rho)
qh<-sig*(1-rho*rho)
for(j in 2:T){
h[j]~dnorm(mu2[j],sig)
mu2[j]<-mu1+rho*h[j-1]}
mu1~dnorm(0,0.0001)
sig~dgamma(0.0001,0.0001)
rho~dunif(-0.999,0.999)}

# Rコード
library(R2WinBUGS)
svdat=read.table("C:/Program Files/R/R-2.7.2/svdat.csv",sep=",",header=T)
T=nrow(svdat)
y=svdat$y
h=svdat$h
data<-list("T","y","h")

inits<-function(){
list(
mu1=rnorm(1,0,0.1),
sig=rgamma(1,0.1,0.1),
rho=runif(1,-0.999,0.999))}

parameters<-c("mu1","sig","rho")

sv.post<-bugs(data,inits,parameters,
model.file="C:/Program Files/R/R-2.7.2/sv.txt",debug=FALSE,
n.iter=10000,n.burnin=1000,n.chains=3,
bugs.directory = "C:/Program Files/WinBUGS14",
working.directory = "C:/Program Files/R/R-2.7.2/")
sv.post
# 図6.12
plot(sv.post)



