# Finding mode of posterior using Newton method in R

I am attempting to approximate the posterior $$\tilde{\pi_{G}}(z|\theta,Y)$$ which is the Gaussian approximation to the full conditional of $$z$$, and in order to do this I need to find the mode $$z^{*} \equiv z^{*}(\theta, Y).$$

Firstly, giving the definition for parameters/variables, we have $$Y_{t} | \log(\lambda_{t}) = x_{t} \sim^{i.i.d} Poi (\log(\lambda _{t})),$$ where $$x_{t} = \log(\lambda_{t})$$ can represent the pattern using the latent autoregressive process, $$x_{t} = \beta_{0} + \beta_{1} x_{t-1} + \eta_{t},$$ where $$\eta_{t} \sim^{i.i.d} N(0, \sigma_{\eta}^{2}).$$ Then taking $$z_t = x_t - \mu$$, we can obtain a nicer form of, $$z_{t} = \beta_{1}z_{t-1} + \eta_{t},$$ which is the $$z$$ used above. This is $$z \sim MVN (0, \Sigma_{z}),$$ and we see that $$z \sim AR(1)$$. To simplify notation we have above, $$\theta = (\beta_0,\beta_1,\sigma^{2}_{\eta})$$ being the parameters, and $$Y$$ is the given observed data.

In order to calculate the mode, we can apply the newton method, $$z^{n} = z^{n-1} + (\nabla^{2}\log(\pi(z^{n-1}|\theta, Y)))^{-1} \nabla \log(\pi(z^{n-1}|\theta,Y)).$$

I have calculate these values as, $$\nabla_{z}\log(\pi(z|\theta,Y)) = \nabla_{z}\log(\pi(z|\theta)) + \nabla_{z}\log(\pi(Y|z)), \quad \text{and}$$ $$\nabla_{z}^{2}\log(\pi(z|\theta,Y)) = \nabla_{z}^{2}\log(\pi(z|\theta)) + \nabla_{z}^{2}\log(\pi(Y|z)),$$ which are simply the gradient and hessian of prior and likelihood function (as we can eliminate any constant terms not involving z, also taking log as it simplifies calculations without changing a stationary point). We can find these through simple calculations as, $$z^{n} = z^{n-1} + (-\Sigma^{-1} - \exp{(\mu)}diag(\exp(z))^{-1} (-\Sigma^{-1}z - \exp{(z)}\exp{(\mu)}+ Y).$$

I have created a code with some synthetic data in order to compute this, however, it seems that with more iterations the magnitude of mode increases, clearly incorrect. Can anyone see what is going wrong? Or any advice?

My code for this is:

# Synthetic data generation
beta0 <- 0
beta1 <- 0.5
varwn <- 1.5
theta<- c(beta0,beta1,varwn)
n = 100
x <- NULL
x[1] <- rnorm(n=1, mean = (beta0/(1-beta1)), sd = sqrt(varwn/(1-beta1^2)))
for (i in 2:n) {
e <- rnorm(n=1, mean = 0, sd = sqrt(varwn))
x[i] <- beta0 +beta1*x[i-1] + e
}
lambda <- exp(x)
Y=rpois(n,lambda = lambda)
data <- Y

# Initialisation of parameters
n     <- length(data)
beta0 <- theta[1]
beta1 <- theta[2]
sigmasqe <- exp(theta[3])
mu       <- (beta0)/(1-beta1)
sigmasq  <- (sigmasqe)/(1-beta1^2)
z<- Matrix(rep(0,n))
znew <- z

# Precision matrix Q = Sigma^-1, (inverse of covariance matrix)
invsigMat            <- Matrix(diag(1 + beta1^2,n))
diag(invsigMat[-1,]) <- -beta1
diag(invsigMat[,-1]) <- -beta1
invsigMat[1,1]       <- 1; invsigMat[n,n] <- 1;
inverseSigMat         = Matrix((1/(sigmasqe))*invsigMat)

#Newton-Raphson, with M iterations
M <- 20
for (i in 1:M){
z = znew