Parameter estimation using the Optim function

Question

I’m having a problem estimating parameters when using the code below:

N=1000
m=matrix(ncol=5,nrow=N)

mu1=1.2
mu2=1.5
phi1=1
phi2=2
rho=0.5
n=100
truevalue=c(mu1,mu2,phi1,phi2,rho)

for (i in 1:N){
U1=rnorm(n,0,1)
U2=rnorm(n,0,1)
z1=U1*((sqrt(1+rho)+sqrt(1-rho))/2)+U2*((sqrt(1+rho)-sqrt(1-rho))/2)
z2=U1*((sqrt(1+rho)-sqrt(1-rho))/2)+U2*((sqrt(1+rho)+sqrt(1-rho))/2)
T1=((phi1*mu1)/(phi1+1))*((1/2)*sqrt(2/phi1)*z1+sqrt(((1/2)*sqrt(2/phi1)*z1)^2+1))^2
T2=((phi2*mu1)/(phi2+1))*((1/2)*sqrt(2/phi2)*z1+sqrt(((1/2)*sqrt(2/phi2)*z1)^2+1))^2
y=T1*(1*(T2>1))

Loglik<-function(theta){

lmu1 <- theta[1]
lmu2 <- theta[2]
lphi1<- theta[3] 
lphi2<- theta[4]
arho <- theta[5]
# Distribuição de Y
gy<-((1/(2*sqrt(2*pi)))*exp((-exp(lphi1)/4)*(sqrt((y*(exp(lphi1)+1))/(exp(lphi1)*exp(lmu1)))-sqrt((exp(lphi1)*exp(lmu1))/(y*(exp(lphi1)+1))))^(2))*(sqrt(exp(lphi1)/2))*(((sqrt(exp(lphi1)+1))/(sqrt(exp(lphi1)*exp(lmu1)*y)))+((sqrt(exp(lphi1)*exp(lmu1)))/(sqrt((exp(lphi1)+1)*y^(3)))))*pnorm(((sqrt(exp(lphi2)*(exp(lphi2)+1))*(exp(lphi2)*exp(lmu2)-(exp(lphi2)+1)))/(sqrt(2)*(exp(lphi2)+1)*sqrt(exp(lphi2)*exp(lmu2)*(1-sin(arho)^2))))+(((sin(arho)*sqrt(exp(lphi1)))/(sqrt(2*(1-sin(arho)^2))))*(sqrt((y*(exp(lphi1)+1))/(exp(lphi1)*exp(lmu1)))-sqrt((exp(lphi1)*exp(lmu1))/(y*(exp(lphi1)+1)))))))
w <- sqrt(exp(lphi2)/2)*(sqrt((exp(lphi2)+1)/(exp(lphi2)*exp(lmu2)))-sqrt(exp(lphi2)*exp(lmu2)/(exp(lphi2)+1)))
fda<- pnorm(w)
g=ifelse((y>0),gy,1)
lv<-sum((y>0)*log(g))+sum((y==0)*log(fda))
return(-lv)
}

grad<-function(theta){
  lmu1 <- theta[1]
  lmu2 <- theta[2]
  lphi1<- theta[3]
  lphi2<- theta[4]
  arho <- theta[5]
  f1<-ifelse((y>0),exp((-(1/4)*exp(lphi1))*((sqrt(y*(exp(lphi1)+1)/(exp(lphi1)*exp(lmu1)))-sqrt(exp(lphi1)*exp(lmu1)/(y*(exp(lphi1)+1))))^2))/(2*sqrt(2*pi)),0)
  f2<-ifelse((y>0),sqrt((1/2)*exp(lphi1)),0)
  f3<-ifelse((y>0),sqrt(exp(lphi1)+1)/sqrt(exp(lphi1)*exp(lmu1)*y)+sqrt(exp(lphi1)*exp(lmu1))/sqrt((exp(lphi1)+1)*y^3),0)
  f4<-ifelse((y>0),pnorm(sqrt(exp(lphi2)*(exp(lphi2)+1))*(exp(lphi2)*exp(lmu2)-exp(lphi2)-1)/(sqrt(2)*(exp(lphi2)+1)*sqrt(exp(lphi2)*exp(lmu2)*(1-sin(arho)^2)))+sqrt(exp(lphi1))*sin(arho)*(sqrt(y*(exp(lphi1)+1)/(exp(lphi1)*exp(lmu1)))-sqrt(exp(lphi1)*exp(lmu1)/(y*(exp(lphi1)+1))))/sqrt(2*(1-sin(arho)^2))),0)
  #Derivada em relação a exp(lmu1)
  df1m1<-ifelse((y>0),-(1/8)*exp(lphi1)*(sqrt(y*(exp(lphi1)+1)/(exp(lphi1)*exp(lmu1)))-sqrt(exp(lphi1)*exp(lmu1)/(y*(exp(lphi1)+1))))*(-(1/2)*y*(exp(lphi1)+1)/(sqrt(y*(exp(lphi1)+1)/(exp(lphi1)*exp(lmu1)))*exp(lphi1)*exp(lmu1)^2)-(1/2)*exp(lphi1)/(sqrt(exp(lphi1)*exp(lmu1)/(y*(exp(lphi1)+1)))*y*(exp(lphi1)+1)))*exp(-(1/4)*exp(lphi1)*(sqrt(y*(exp(lphi1)+1)/(exp(lphi1)*exp(lmu1)))-sqrt(exp(lphi1)*exp(lmu1)/(y*(exp(lphi1)+1))))^2)*sqrt(2)/sqrt(pi),0)
  df2m1<-0
  df3m1<-ifelse((y>0),-(1/2)*sqrt(exp(lphi1)+1)*exp(lphi1)*y/(exp(lphi1)*exp(lmu1)*y)^(3/2)+(1/2)*exp(lphi1)/(sqrt(exp(lphi1)*exp(lmu1))*sqrt((exp(lphi1)+1)*y^3)),0)
  df4m1<-ifelse((y>0),((1/2)*exp(-(1/2)*((1/2)*sqrt(exp(lphi2)*(exp(lphi2)+1))*(exp(lphi2)*exp(lmu2)-exp(lphi2)-1)*sqrt(2)/sqrt(exp(lphi2)*exp(lmu2)*(1-sin(arho)^2))+sqrt(exp(lphi1))*sin(arho)*(sqrt(y*(exp(lphi1)+1)/(exp(lphi1)*exp(lmu1)))-sqrt(exp(lphi1)*exp(lmu1)/(y*(exp(lphi1)+1))))/sqrt(2-2*sin(arho)^2))^2)*sqrt(2)/sqrt(pi))*(sqrt(exp(lphi1))*sin(arho)*(-(1/2)*y*(exp(lphi1)+1)/(sqrt(y*(exp(lphi1)+1)/(exp(lphi1)*exp(lmu1)))*exp(lphi1)*exp(lmu1)^2)-(1/2)*exp(lphi1)/(sqrt(exp(lphi1)*exp(lmu1)/(y*(exp(lphi1)+1)))*y*(exp(lphi1)+1)))/sqrt(2-2*sin(arho)^2)),0)
  #Derivada em relação a exp(lmu2)
  df1m2<-0
  df2m2<-0
  df3m2<-0
  df4m2<-ifelse((y>0),((1/2)*exp(-(1/2)*((1/2)*sqrt(exp(lphi2)*(exp(lphi2)+1))*(exp(lphi2)*exp(lmu2)-exp(lphi2)-1)*sqrt(2)/sqrt(exp(lphi2)*exp(lmu2)*(1-sin(arho)^2))+sqrt(exp(lphi1))*sin(arho)*(sqrt(y*(exp(lphi1)+1)/(exp(lphi1)*exp(lmu1)))-sqrt(exp(lphi1)*exp(lmu1)/(y*(exp(lphi1)+1))))/sqrt(2-2*sin(arho)^2))^2)*sqrt(2)/sqrt(pi))*((1/2)*sqrt(exp(lphi2)*(exp(lphi2)+1))*exp(lphi2)*sqrt(2)/sqrt(exp(lphi2)*exp(lmu2)*(1-sin(arho)^2))-(1/4)*sqrt(exp(lphi2)*(exp(lphi2)+1))*(exp(lphi2)*exp(lmu2)-exp(lphi2)-1)*sqrt(2)*exp(lphi2)*(1-sin(arho)^2)/(exp(lphi2)*exp(lmu2)*(1-sin(arho)^2))^(3/2)),0)
  #Derivada em relação a exp(lphi1)
  df1ph1<-ifelse((y>0),(1/4)*(-(1/4)*(sqrt(y*(exp(lphi1)+1)/(exp(lphi1)*exp(lmu1)))-sqrt(exp(lphi1)*exp(lmu1)/(y*(exp(lphi1)+1))))^2-(1/2)*exp(lphi1)*(sqrt(y*(exp(lphi1)+1)/(exp(lphi1)*exp(lmu1)))-sqrt(exp(lphi1)*exp(lmu1)/(y*(exp(lphi1)+1))))*((1/2)*(y/(exp(lphi1)*exp(lmu1))-y*(exp(lphi1)+1)/(exp(lphi1)^2*exp(lmu1)))/sqrt(y*(exp(lphi1)+1)/(exp(lphi1)*exp(lmu1)))-(1/2)*(exp(lmu1)/(y*(exp(lphi1)+1))-exp(lphi1)*exp(lmu1)/(y*(exp(lphi1)+1)^2))/sqrt(exp(lphi1)*exp(lmu1)/(y*(exp(lphi1)+1)))))*exp(-(1/4)*exp(lphi1)*(sqrt(y*(exp(lphi1)+1)/(exp(lphi1)*exp(lmu1)))-sqrt(exp(lphi1)*exp(lmu1)/(y*(exp(lphi1)+1))))^2)*sqrt(2)/sqrt(pi),0)
  df2ph1<-ifelse((y>0),(1/4)*sqrt(2)/sqrt(exp(lphi1)),0)
  df3ph1<-ifelse((y>0),1/(2*sqrt(exp(lphi1)+1)*sqrt(exp(lphi1)*exp(lmu1)*y))-(1/2)*sqrt(exp(lphi1)+1)*exp(lmu1)*y/(exp(lphi1)*exp(lmu1)*y)^(3/2)+(1/2)*exp(lmu1)/(sqrt(exp(lphi1)*exp(lmu1))*sqrt((exp(lphi1)+1)*y^3))-(1/2)*sqrt(exp(lphi1)*exp(lmu1))*y^3/((exp(lphi1)+1)*y^3)^(3/2),0)
  df4ph1<-ifelse((y>0),((1/2)*exp(-(1/2)*((1/2)*sqrt(exp(lphi2)*(exp(lphi2)+1))*(exp(lphi2)*exp(lmu2)-exp(lphi2)-1)*sqrt(2)/sqrt(exp(lphi2)*exp(lmu2)*(1-sin(arho)^2))+sqrt(exp(lphi1))*sin(arho)*(sqrt(y*(exp(lphi1)+1)/(exp(lphi1)*exp(lmu1)))-sqrt(exp(lphi1)*exp(lmu1)/(y*(exp(lphi1)+1))))/sqrt(2-2*sin(arho)^2))^2)*sqrt(2)/sqrt(pi))*((1/2)*sin(arho)*(sqrt(y*(exp(lphi1)+1)/(exp(lphi1)*exp(lmu1)))-sqrt(exp(lphi1)*exp(lmu1)/(y*(exp(lphi1)+1))))/(sqrt(exp(lphi1))*sqrt(2-2*sin(arho)^2))+sqrt(exp(lphi1))*sin(arho)*((1/2)*(y/(exp(lphi1)*exp(lmu1))-y*(exp(lphi1)+1)/(exp(lphi1)^2*exp(lmu1)))/sqrt(y*(exp(lphi1)+1)/(exp(lphi1)*exp(lmu1)))-(1/2)*(exp(lmu1)/(y*(exp(lphi1)+1))-exp(lphi1)*exp(lmu1)/(y*(exp(lphi1)+1)^2))/sqrt(exp(lphi1)*exp(lmu1)/(y*(exp(lphi1)+1))))/sqrt(2-2*sin(arho)^2)),0)
  #Derivada em relação a exp(lphi2)
  df1ph2<-0
  df2ph2<-0
  df3ph2<-0
  df4ph2<-ifelse((y>0),((1/2)*exp(-(1/2)*((1/2)*sqrt(exp(lphi2)*(exp(lphi2)+1))*(exp(lphi2)*exp(lmu2)-exp(lphi2)-1)*sqrt(2)/sqrt(exp(lphi2)*exp(lmu2)*(1-sin(arho)^2))+sqrt(exp(lphi1))*sin(arho)*(sqrt(y*(exp(lphi1)+1)/(exp(lphi1)*exp(lmu1)))-sqrt(exp(lphi1)*exp(lmu1)/(y*(exp(lphi1)+1))))/sqrt(2-2*sin(arho)^2))^2)*sqrt(2)/sqrt(pi))*((1/2)*sin(arho)*(sqrt(y*(exp(lphi1)+1)/(exp(lphi1)*exp(lmu1)))-sqrt(exp(lphi1)*exp(lmu1)/(y*(exp(lphi1)+1))))/(sqrt(exp(lphi1))*sqrt(2-2*sin(arho)^2))+sqrt(exp(lphi1))*sin(arho)*((1/2)*(y/(exp(lphi1)*exp(lmu1))-y*(exp(lphi1)+1)/(exp(lphi1)^2*exp(lmu1)))/sqrt(y*(exp(lphi1)+1)/(exp(lphi1)*exp(lmu1)))-(1/2)*(exp(lmu1)/(y*(exp(lphi1)+1))-exp(lphi1)*exp(lmu1)/(y*(exp(lphi1)+1)^2))/sqrt(exp(lphi1)*exp(lmu1)/(y*(exp(lphi1)+1))))/sqrt(2-2*sin(arho)^2)),0)
  #Derivada em relação a sin(arho)
  df1rh<-0
  df2rh<-0
  df3rh<-0
  df4rh<-ifelse((y>0),((1/2)*exp(-(1/2)*((1/2)*sqrt(exp(lphi2)*(exp(lphi2)+1))*(exp(lphi2)*exp(lmu2)-exp(lphi2)-1)*sqrt(2)/sqrt(exp(lphi2)*exp(lmu2)*(1-sin(arho)^2))+sqrt(exp(lphi1))*sin(arho)*(sqrt(y*(exp(lphi1)+1)/(exp(lphi1)*exp(lmu1)))-sqrt(exp(lphi1)*exp(lmu1)/(y*(exp(lphi1)+1))))/sqrt(2-2*sin(arho)^2))^2)*sqrt(2)/sqrt(pi))*((1/2)*sqrt(exp(lphi2)*(exp(lphi2)+1))*(exp(lphi2)*exp(lmu2)-exp(lphi2)-1)*sqrt(2)*exp(lphi2)*exp(lmu2)*sin(arho)/(exp(lphi2)*exp(lmu2)*(1-sin(arho)^2))^(3/2)+sqrt(exp(lphi1))*(sqrt(y*(exp(lphi1)+1)/(exp(lphi1)*exp(lmu1)))-sqrt(exp(lphi1)*exp(lmu1)/(y*(exp(lphi1)+1))))/sqrt(2-2*sin(arho)^2)+2*sqrt(exp(lphi1))*sin(arho)^2*(sqrt(y*(exp(lphi1)+1)/(exp(lphi1)*exp(lmu1)))-sqrt(exp(lphi1)*exp(lmu1)/(y*(exp(lphi1)+1))))/(2-2*sin(arho)^2)^(3/2)),0)
  grad1<-sum(df1m1*f2*f3*f4+f1*f2*df3m1*f4+f1*f2*f3*df4m1)
  grad2<-sum(f1*f2*f3*df4m2)
  grad3<-sum(df1ph1*f2*f3*f4+f1*df2ph1*f3*f4+f1*f2*df3ph1*f4+f1*f2*f3*df4ph1)
  grad4<-sum(f1*f2*f3*df4ph2)
  grad5<-sum(f1*f2*f3*df4rh)
  c(-grad1,-grad2,-grad3,-grad4,-grad5)
}

lmu1_0=log(2)
lmu2_0=log(1.5)
lphi1_0=log(3)
lphi2_0=log(3)
arho_0=asin(0.6)
start=c(lmu1_0,lmu2_0,lphi1_0,lphi2_0,arho_0)

op=optim(par=start, Loglik, grad, method="BFGS", hessian=F, control= list(maxit=1000000))
    m[i,]=c(exp(op$par[1]),exp(op$par[2]),exp(op$par[3]),exp(op$par[4]),sin(op$par[5    ]))
}

Does anyone have any suggestion to improve the estimation of the parameters, also imprementei the gradient and used the same to run the function Optim. To facilitate the analysis a reparamerization of the parameters was made using functions 1 to 1 (logarithm and sine). The parameters mu1, mu2,phi1 and phi2 must be positive numbers and rhois a number between -1 and 1. I appreciate the help!!!

Answer 1

It was not clear to me what it means to "improve convergence". I ran your code on my machine and the function converged.

If "improving convergence" means obtaining similar results in fewer interactions, I suggest using the nlm function. Her problem is being more sensitive to the shape of the likelihood function.

If "improving convergence" means getting values closer to the values used in the simulation, increase the number of iterations beyond 1e6 or use the abstol argument with a very small value.

However, I see no reason to do this. The estimates obtained using your command seem to be excellent.