##################################
# AR(1)
##################################
set.seed(12345)
n = 100
alpha = 0.1
phit=0.9
sig2 = 0.1
y = rep(0,2*n)
for (t in 2:(2*n))
  y[t]=rnorm(1,alpha+phit*y[t-1],sqrt(sig2))
y = y[(n+1):(2*n)]

pdf(file="ar1-graph1.pdf",width=8,height=6)
ts.plot(y) 
dev.off()

b = 0
B = 1

N1 = 10000
h = phis1[2]-phis1[1]
phis1 = seq(-1,1,length=N1)
phis2 = seq(-3,3,by=h)
N2 = length(phis2)
like1 = rep(0,N1)
like2 = rep(0,N2)
like1n = rep(0,N1)
like2n = rep(0,N2)

for (i in 1:N1){
  like1[i] = dnorm(y[1],alpha/(1-phis1[i]),sqrt(sig2/(1-phis1[i]^2)))
  like1n[i] = prod(dnorm(y[2:n],phis1[i]*y[1:(n-1)],sqrt(sig2)))
}
for (i in 1:N2){
  like2[i] = dnorm(y[1],alpha+phis2[i]*b,sqrt(phis2[i]^2*B+sig2))
  like2n[i] = prod(dnorm(y[2:n],phis2[i]*y[1:(n-1)],sqrt(sig2)))
}
like1 = like1/sum(like1)/h
like2 = like2/sum(like2)/h
like1n = like1n/sum(like1n)/h
like2n = like2n/sum(like2n)/h
post1 = like1*like1n
post1 = post1/sum(post1)/h
post2 = like2*like2n
post2 = post2/sum(post2)/h

pdf(file="ar1-graph2.pdf",width=11,height=7)
par(mfrow=c(1,2))
plot(phis1,like1,type="l",xlab=expression(phi),ylab="Likelihood",ylim=range(like1,like2))
title(substitute(paste(y[1],"|",phi," ~ N[",alpha/(1-phi),",",sigma^2/(1-phi^2),"]")))
legend("topleft",legend=substitute(paste("True: ",alpha,"=",a," , ",phi,"=",b," , ",sigma^2," = ",s2," , ",y[1]," = ",m),list(a=alpha,b=phit,s2=sig2,m=round(y[1],3))))

plot(phis2,like2,type="l",xlab=expression(phi),ylab="Likelihood",ylim=range(like1,like2))
title(substitute(paste(y[1],"|",phi," ~ N[",alpha+phi * b,",",phi^2 * B + sigma^2,"]")))
legend("topleft",legend=paste("Hyperparameters: b=",b,", B=",B,sep=""))
dev.off()




pdf(file="ar1-graph3.pdf",width=11,height=7)

par(mfrow=c(1,2))
plot(phis1,like1,type="l",xlab=expression(phi),ylab="Likelihood",ylim=range(like1,like2,like1n,like2n))
title(substitute(paste(y[1],"|",phi," ~ N[",alpha/(1-phi),",",sigma^2/(1-phi^2),"]")))
lines(phis1,like1n,col=2)
legend("topleft",legend=c("Likelihood for 1 obs.",paste("Likelihood for ",n-1," obs.",sep="")),col=1:2,lty=1)


plot(phis2,like2,type="l",xlab=expression(phi),ylab="Likelihood",ylim=range(like1,like2,like1n,like2n))
title(substitute(paste(y[1],"|",phi," ~ N[",alpha+phi * b,",",phi^2 * B + sigma^2,"]")))
lines(phis2,like2n,col=2)
dev.off()


pdf(file="ar1-graph4.pdf",width=8,height=6)
par(mfrow=c(1,1))
plot(phis1,post1,type="l",xlab=expression(phi),ylab="Likelihood",ylim=range(post1,post2),xlim=c(0.85,1.05))
lines(phis2,post2,col=2)
abline(v=1,lty=3)
legend("topleft",legend=c("Stationarity constrained","Unconstrained"),col=1:2,lty=1)
title("Posterior distributions")
dev.off()




##################################
# AR(2)
##################################
set.seed(21356)
n = 100
alpha = 0.1
phi1t=0.25
phi2t=0.25
sig2 = 0.1
y = rep(0,2*n)
for (t in 3:(2*n))
  y[t]=rnorm(1,alpha+phi1t*y[t-1]+phi2t*y[t-2],sqrt(sig2))
y = y[(n+1):(2*n)]

ts.plot(y) 
 
b = c(0.2,0.2)
B = diag(1000,2)

mm =rep(0,2)
VV  = diag(0,2)
d     = rep(0,2)
D    = diag(0,2)

install.packages("mvtnorm")
library("mvtnorm")

N = 100
phi1 = seq(0,0.5,length=N)
phi2 = seq(0,0.5,length=N)
h = phi1[2]-phi1[1]
like1 = matrix(0,N,N)
like2 = matrix(0,N,N)
like3 = matrix(0,N,N)
for (i in 1:N)
for (j in 1:N){
m       = rep(alpha/(1-phi1[i]-phi2[j]),2)
V[1,2] = phi1[i]/((1+phi2[j])*((1-phi2[j])^2-phi1[i]^2))*sig2
V[2,1] = V[1,2]
V[1,1] = (1-phi2[i])/((1+phi2[j])*((1-phi2[j])^2-phi1[i]^2))*sig2
V[2,2] =V [1,1]

d[1]     = alpha+phi1[i]*b[1]+phi2[j]*b[2]
d[2]     = alpha*(1+phi1[i]) + b[1]*(phi1[i]^2+phi2[j])+b[2]*phi1[i]*phi2[j]
D[1,1] = phi1[i]^2*B[1,1]+phi2[j]^2*B[2,2]+2*phi1[i]*phi2[j]*B[1,2]+sig2
D[2,2] = phi1[i]^2*D[1,1]+phi2[j]*B[1,1]+2*phi1[i]*phi2[j]*(phi1[i]*B[1,1]+phi2[j]*B[1,2])+sig2
D[1,2] = phi1[i]*phi2[j]*B[1,1]+phi1[i]^2*(phi1[i]*B[1,1]+phi2[j]*B[1,2])+
           phi1[i]*phi2[j]*(phi1[i]*B[1,2]+phi2[j]*B[2,2])+phi2[j]^2*B[1,2]
D[2,1] = D[1,2]
like1[i,j]=dmvnorm(c(y[1],y[2]),m,V)
like2[i,j]=dmvnorm(c(y[1],y[2]),d,D)
like3[i,j]=prod(dnorm(y[3:n],alpha+phi1[i]*y[2:(n-1)]+phi2[j]*y[1:(n-2)],sqrt(sig2)))
}
like1 = like1/sum(like1)/(h^2)
like2 = like2/sum(like2)/(h^2)
like3 = like3/sum(like3)/(h^2)

par(mfrow=c(1,3))
contour(phi1,phi2,like1,xlab=expression(phi[1]),ylab=expression(phi[2]))
points(phi1t,phi2t,col=2,pch=16)
title("Likelihood based on y1,y2\nStationarity imposed")
contour(phi1,phi2,like2,xlab=expression(phi[1]),ylab=expression(phi[2]))
points(phi1t,phi2t,col=2,pch=16)
title("Likelihood based on y1,y2\nStationarity not imposed")
contour(phi1,phi2,like3,xlab=expression(phi[1]),ylab=expression(phi[2]))
points(phi1t,phi2t,col=2,pch=16)
title("Likelihood based on y3,...,y100")

par(mfrow=c(1,3))
contour(phi1,phi2,like3,xlab=expression(phi[1]),ylab=expression(phi[2]),drawlabels=FALSE)
points(phi1t,phi2t,col=2,pch=16)
contour(phi1,phi2,like1*like3,xlab=expression(phi[1]),ylab=expression(phi[2]),drawlabels=FALSE)
points(phi1t,phi2t,col=2,pch=16)
contour(phi1,phi2,like2*like3,xlab=expression(phi[1]),ylab=expression(phi[2]),drawlabels=FALSE)
points(phi1t,phi2t,col=2,pch=16)