estou tendo um problema na estimação de parâmetros ao usar o código abaixo:
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 ]))
}
Será que alguém tem alguma sugestão para melhorar a estimação dos parâmetros, imprementei também o gradiente e utilizei o mesmo para rodar a função optim. Para facilitar a análise foi feito uma reparamerização dos parâmetros usando funções 1 a 1 (logaritmo e seno). Os parâmetros mu1
, mu2
,phi1
e phi2
devem ser números positivos e rho
é um número entre -1 e 1. Desde já agradeço a ajuda!!!