1

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!!!

1 Resposta 1

1

Não ficou claro, para mim, o que significa "melhorar a convergência". Rodei o teu código na minha máquina e a função optim convergiu.

Se "melhorar a convergência" significar obter resultados similares em menos interações, sugiro utilizar a função nlm. O problema dela é ser mais sensível à forma da função de verossimilhança.

Se "melhorar a convergência" significar obter valores mais próximos dos valores utilizados na simulação, aumente o número de iterações para além de 1e6 ou utilize o argumento abstol com um valor bem pequeno.

Entretanto, não vejo motivo para fazer isto. As estimativas obtidas utilizando o teu comando parecem ser excelentes.

4
  • Olá amigo, por "melhorar a convergência" quero dizer obter valores mais próximos dos valores utilizados na simulação sim. Ao rodar o código meu rho tem dado quase sempre igual a 1, quando na simulação usei 0.5, além disso, os demais parâmetros tenho obtido com um viés alto. Ou seja, a estimação não está excelente, talvez não esteja incorreta, mas gostaria de ver alguma forma de melhorar!
    – fsbmat
    Commented 5/07/2016 às 17:00
  • Já tentei aumentar o número de iterações e usar o argumento abstol com um valor pequeno, mas não melhorou a estimação!
    – fsbmat
    Commented 5/07/2016 às 17:11
  • 1
    Para quem viu esta dúvida, imprementei também o gradiente, mas o problema de estimação persiste, agora parece que a função optim não está considerando os valores iniciais ou não está resubstituindo os valores estimados em cada passo, se alguém tiver alguma sugestão fico muito grato, alterei o primeiro código postado agora o mesmo segue com o gradiente!
    – fsbmat
    Commented 6/07/2016 às 22:02
  • o problema do meu código era o gradiente, o mesmo não estava correto. Continuo com alguns outros problemas, mas em relação ao erro postado aqui o problema era o gradiente! Agradeço a ajuda @grandeabobora, obrigado! O código acima está com o gradiente corrigido!!!
    – fsbmat
    Commented 7/07/2016 às 13:07

Você deve fazer log-in para responder a esta pergunta.

Esta não é a resposta que você está procurando? Pesquise outras perguntas com a tag .