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)