Parameter estimation using the optimum function

Question

Parameter estimation using the optimum function

Navigation

#1 by (1 votes)

1

I'm having a problem in parameter estimation 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 suggestions for improving the parameter estimation, I have also imprinted the gradient and used it to run the optimum function. To facilitate the analysis a reparamerization of the parameters was done using functions 1 to 1 (log and sine). The parameters mu1 , mu2 , phi1 and phi2 must be positive numbers and rho is a number between -1 and 1. From now on thanks for help !!!

r

asked by anonymous 04.07.2016 / 23:44

1 answer

How to configure JBoss AS 7.1 with Demoiselle 2.4.1 to operate with distributed transactions? How to configure the remote GIT in Amazon AWS?

score 1 · Answer 1

It has not been clear to me what "improving convergence" means. I circled your code on my machine and the optimum function converged.

If "improving convergence" means achieving similar results in less interactions, I suggest using the nlm function. Her problem is to be more sensitive to the form 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 appear to be excellent.