set.seed(234345)
n = 200
sig = 1
phi = c(1.5,-0.75)
y = rep(0,2*n)
for (t in 3:(2*n))
  y[t] = rnorm(1,sum(phi*c(y[t-1],y[t-2])),sig)
y = y[(n+1):(2*n)]


pdf(file="AR2-bayes.pdf",width=10,height=8)
par(mfrow=c(1,1))
ts.plot(y)

par(mfrow=c(1,2))
acf(y)
pacf(y)

# MLE/OLS 

mleARp = function(y,n,p){
  n1 = n-p
  y1 = y[(p+1):n]
  X = matrix(0,n1,p)
  for (i in 1:p)
    X[,i] = y[(p-i+1):(n-i)]
  iXtX = solve(t(X)%*%X)
  phihat = iXtX%*%t(X)%*%y1
  yhat = X%*%phihat
  sig2hat = sum((y1-yhat)^2)/(n1-p)
  sephihat = sqrt(sig2hat*diag(iXtX))
  ttests = phihat/sephihat 
  pvalues = 2*pt(-abs(ttests),n1-p)
  regtable = cbind(round(phihat,5),round(sephihat,5),round(ttests,3),round(pvalues,4))
  return(regtable)
}


bayesARp = function(y,n,p,b0,B0,n0,s02){
  n1 = n-p
  y1 = y[(p+1):n]
  X = matrix(0,n1,p)
  for (i in 1:p)
    X[,i] = y[(p-i+1):(n-i)]
  B1 = solve(B0+t(X)%*%X)
  b1 = B1%*%(solve(B0)%*%b0+t(X)%*%y1)
  n1 = n-2+n0
  s12 = (n0*s02+t(y1-X%*%b1)%*%y1+ t(b0-b1)%*%solve(B0)%*%b0)/n1
  sephi = sqrt(s12*diag(B1))
  ttests = b1/sephi
  pvalues = 2*pt(-abs(ttests),n1)
  regtable = cbind(round(b1,5),round(sephi,5),round(ttests,3),round(pvalues,4))
  return(list(b=b1,B=B1,n=n1,s2=s12,table=regtable))
}


mleARp(y,n,2)

b0 = rep(0,2)
B0 = diag(1,2)
#b0 = c(-3,3)
#B0 = diag(0.01,2)
n0 = 5
s02 = 1
runbayes = bayesARp(y,n,2,b0,B0,n0,s02)
runbayes$table

n1 = runbayes$n
s12 = runbayes$s2
b1 = runbayes$b
B1 = runbayes$B


library(MASS)
set.seed(2325)
M = 10000
draw.sig20 = 1/rgamma(M,n0/2,n0*s02/2)
draw.sig2 = 1/rgamma(M,n1/2,n1*s12/2)
draw.phi0 = matrix(0,M,2)
draw.phi = matrix(0,M,2)
for (i in 1:M){
  draw.phi0[i,] = mvrnorm(1,b0,sqrt(draw.sig20[i])*B0)
  draw.phi[i,] = mvrnorm(1,b1,sqrt(draw.sig2[i])*B1)
}
draw.sig0 = sqrt(draw.sig20)
draw.sig = sqrt(draw.sig2)


par(mfrow=c(1,1))
plot(draw.phi0[,1],draw.phi0[,2],xlab=expression(phi[1]),ylab=expression(phi[2]),
pch=16,col=grey(0.5),xlim=c(-4,4),ylim=c(-4,4))
points(draw.phi[,1],draw.phi[,2],pch=16)
points(phi[1],phi[2],col=2,pch=16)
segments(0,1,2,-1,col=4)
segments(0,1,-2,-1,col=4)
segments(-2,-1,2,-1,col=4)
abline(h=0,lty=2)
abline(v=0,lty=2)


par(mfrow=c(2,3))
plot(draw.phi0[,1],draw.phi0[,2],xlab=expression(phi[1]),ylab=expression(phi[2]),pch=16,col=grey(0.5))
points(draw.phi[,1],draw.phi[,2],pch=16)
points(phi[1],phi[2],col=2,pch=16)
segments(0,1,2,-1,col=4)
segments(0,1,-2,-1,col=4)
segments(-2,-1,2,-1,col=4)
abline(h=0,lty=2)
abline(v=0,lty=2)

plot(draw.phi0[,1],draw.sig0,xlab=expression(phi[1]),ylab=expression(sigma),pch=16,col=grey(0.5))
points(draw.phi[,1],draw.sig,pch=16)
points(phi[1],sig,col=2,pch=16)

plot(draw.phi0[,2],draw.sig0,xlab=expression(phi[2]),ylab=expression(sigma),pch=16,col=grey(0.5))
points(draw.phi[,2],draw.sig,pch=16)
points(phi[2],sig,col=2,pch=16)


plot(density(draw.phi0[,1]),xlim=range(draw.phi[,1],draw.phi0[,1]),main=expression(phi[1]),xlab="")
lines(density(draw.phi[,1]),col=grey(0.5))
abline(v=phi[1],col=2)

plot(density(draw.phi0[,2]),xlim=range(draw.phi[,2],draw.phi0[,2]),main=expression(phi[2]),xlab="")
lines(density(draw.phi[,2]),col=grey(0.5))
abline(v=phi[2],col=2)


plot(density(draw.sig0),xlim=range(draw.sig,draw.sig0),main=expression(sigma),xlab="")
lines(density(draw.sig),col=grey(0.5))
abline(v=sig,col=2)

set.seed(132345)
M = 10000
H = 50
phis = draw.phi
sigs = draw.sig
ypred = matrix(c(y[n],y[n-1]),2,M)
yf = matrix(0,M,H)
for (h in 1:H){
	mean = diag(phis%*%ypred)
	yf[,h] = rnorm(M,mean,sigs)
    ypred = rbind(yf[,h],ypred[1,])
}

M1 = 200
ind = sample(1:M,size=M1,replace=FALSE)
par(mfrow=c(1,1))
plot(y,xlim=c(1,n+H),xlab="time",ylab="",ylim=range(y,yf),type="l")
lines((n+1):(n+H),yf[ind[1],],col=grey(0.5))
for (i in 2:M1)
lines((n+1):(n+H),yf[ind[i],],col=grey(0.5))
qyf = apply(yf,2,quantile,c(0.05,0.5,0.95))
lines((n+1):(n+H),qyf[1,],lwd=2)
lines((n+1):(n+H),qyf[2,],lwd=2)
lines((n+1):(n+H),qyf[3,],lwd=2)
abline(v=n,lty=2)
dev.off()