Erros padrão das estimativas de distribuição hiperbólica usando o método delta?

8

Quero calcular os erros padrão de uma distribuição hiperbólica ajustada.

Na minha notação, a densidade é dada por

H(l;α,β,μ,δ)=α2β22αδK1(δα2β2)exp(αδ2+(lμ)2+β(lμ))
Estou usando o pacote HyperbolicDistr em R. Estimo os parâmetros através do seguinte comando:
hyperbFit(mydata,hessian=TRUE)

Isso me dá uma parametrização errada. Eu mudo para a minha parametrização desejada com o hyperbChangePars(from=1,to=2,c(mu,delta,pi,zeta))comando Então eu quero ter os erros padrão de minhas estimativas, posso obtê-lo para a parametrização errada com o summarycomando Mas isso me dá os erros padrão para a outra parametrização. De acordo com este tópico , tenho que usar o método delta ( não quero usar autoinicialização ou validação cruzada).

O hyperbFitcódigo está aqui . E o hyperbChangeParsestá aqui . Portanto, eu sei que μ e δ permanecem os mesmos. Portanto, também os erros padrão são os mesmos, certo?

Para transformar e ζ em α e β , preciso da relação entre eles. De acordo com o código, isso é feito da seguinte maneira:πζαβ

alpha <- zeta * sqrt(1 + hyperbPi^2) / delta
beta <- zeta * hyperbPi / delta

Então, como preciso codificar o método delta para obter os erros padrão desejados?

EDIT: eu estou usando esses dados . Primeiro, realizo o método delta de acordo com este segmento.

# fit the distribution

hyperbfitdb<-hyperbFit(mydata,hessian=TRUE)
hyperbChangePars(from=1,to=2,hyperbfitdb$Theta)
summary(hyperbfitdb)

summary(hyperbfitdb) fornece a seguinte saída:

Data:      mydata 
Parameter estimates:
        pi           zeta         delta           mu    
    0.0007014     1.3779503     0.0186331    -0.0001352 
  ( 0.0938886)  ( 0.9795029)  ( 0.0101284)  ( 0.0035774)
Likelihood:         615.992 
Method:             Nelder-Mead 
Convergence code:   0 
Iterations:         315 

e hyperbChangePars(from=1,to=2,hyperbfitdb$Theta)fornece a seguinte saída:

   alpha.zeta     beta.zeta   delta.delta         mu.mu 
73.9516898823  0.0518715378  0.0186331187 -0.0001352342 

agora eu defino as variáveis ​​da seguinte maneira:

pi<-0.0007014 
lzeta<-log(1.3779503)
ldelta<-log(0.0186331)

Agora corro o código (segunda edição) e obtenho o seguinte resultado:

> se.alpha
         [,1]
[1,] 13.18457
> se.beta
        [,1]
[1,] 6.94268

Isso está correto? Gostaria de saber o seguinte: Se eu usar um algoritmo de autoinicialização da seguinte maneira:

B = 1000 # number of bootstraps

alpha<-NA
beta<-NA
delta<-NA
mu<-NA


# Bootstrap
for(i in 1:B){
  print(i)
  subsample = sample(mydata,rep=T)
  hyperboot <- hyperbFit(subsample,hessian=FALSE)
  hyperboottransfparam<- hyperbChangePars(from=1,to=2,hyperboot$Theta)
  alpha[i]    = hyperboottransfparam[1]
  beta[i]    = hyperboottransfparam[2]
  delta[i] = hyperboottransfparam[3]
  mu[i] = hyperboottransfparam[4]

}
# hist(beta,breaks=100,xlim=c(-200,200))
sd(alpha)
sd(beta)
sd(delta)
sd(mu)

Eu recebo 119.6por sd(alpha)e 35.85para sd(beta). Os resultados são muito diferentes? Há algum erro ou qual é o problema aqui?

Jen Bohold
fonte

Respostas:

10

hyperbPiπsummaryhyperbFitVar(X)=SE(X)2gα(ζ,π,δ)gβ(ζ,π,δ)αβ

gα(ζ,π,δ)=ζ1+π2δgβ(ζ,π,δ)=ζπδ
α
ζgα(ζ,π,δ)=1+π2δπgα(ζ,π,δ)=πζ1+π2δδgα(ζ,π,δ)=1+π2ζδ2
β
ζgβ(ζ,π,δ)=πδπgβ(ζ,π,δ)=ζδδgβ(ζ,π,δ)=πζδ2

α

Var(α)1+π2δ2Var(ζ)+π2ζ2(1+π2)δ2Var(π)+(1+π2)ζ2δ4Var(δ)+2×[πζδ2Cov(π,ζ)(1+π2)ζδ3Cov(δ,ζ)πζ2δ3Cov(δ,π)]
β

Var(β)π2δ2Var(ζ)+ζ2δ2Var(π)+π2ζ2δ4Var(δ)+2×[πζδ2Cov(π,ζ)π2ζδ3Cov(δ,ζ)πζ2δ3Cov(π,δ)]

Codificação em R

Dαβζ,π,δΣ3×3ζ,π,δvcov(my.hyperbFit)my.hyperbFitα

Var(α)DαΣDα
β

Em R, isso pode ser facilmente codificado como este:

#-----------------------------------------------------------------------------
# The row vector D of the partial derivatives for alpha
#-----------------------------------------------------------------------------

D.alpha <- matrix(
  c(
    sqrt(1+pi^2)/delta,                 # differentiate wrt zeta
    ((pi*zeta)/(sqrt(1+pi^2)*delta)),   # differentiate wrt pi
    -(sqrt(1+pi^2)*zeta)/(delta^2)      # differentiate wrt delta
  ),
  ncol=3)

#-----------------------------------------------------------------------------
# The row vector D of the partial derivatives for beta
#-----------------------------------------------------------------------------

D.beta <- matrix(
  c(
    (pi/delta),            # differentiate wrt zeta
    (zeta/delta),          # differentiate wrt pi
    -((pi*zeta)/delta^2)   # differentiate wrt delta
  ),
  ncol=3)

#-----------------------------------------------------------------------------
# Calculate the approximations of the variances for alpha and beta
# "sigma" denotes the 3x3 covariance matrix
#-----------------------------------------------------------------------------

var.alpha <- D.alpha %*% sigma %*% t(D.alpha) 
var.beta <- D.beta %*% sigma %*% t(D.beta)

#-----------------------------------------------------------------------------
# The standard errors are the square roots of the variances
#-----------------------------------------------------------------------------

se.alpha <- sqrt(var.alpha)
se.beta <- sqrt(var.beta)

log(ζ)log(δ)

ζ=log(ζ)δ=log(δ)ζδ

gα(ζ,π,δ)=exp(ζ)1+π2exp(ζ)gβ(ζ,π,δ)=exp(ζ)πexp(δ)
α
ζgα(ζ,π,δ)=1+π2exp(δ+ζ)πgα(ζ,π,δ)=πexp(δ+ζ)1+π2δgα(ζ,π,δ)=1+π2exp(δ+ζ)
β
ζgβ(ζ,π,δ)=πexp(δ+ζ)πgβ(ζ,π,δ)=exp(δ+ζ)δgβ(ζ,π,δ)=πexp(δ+ζ)
α
Var(α)(1+π2)exp(2δ+2ζ)Var(ζ)+π2exp(2δ+2ζ)1+π2Var(π)+(1+π2)exp(2δ+2ζ)Var(δ)+2×[πexp(2δ+2ζ)Cov(π,ζ)(1+π2)exp(2δ+2ζ)Cov(δ,ζ)πexp(2δ+2ζ)Cov(δ,π)]
β
Var(β)π2exp(2δ+2ζ)Var(ζ)+exp(2δ+2ζ)Var(π)+π2exp(2δ+2ζ)Var(δ)+2×[πexp(2δ+2ζ)Cov(π,ζ)π2exp(2δ+2ζ)Cov(δ,ζ)πexp(2δ+2ζ)Cov(δ,π)]

Codificação em R2

sigmaζ=log(ζ)δ=log(δ)ζδ

#-----------------------------------------------------------------------------
# The row vector D of the partial derivatives for alpha
#-----------------------------------------------------------------------------

D.alpha <- matrix(
  c(
    sqrt(1+pi^2)*exp(-ldelta + lzeta),            # differentiate wrt lzeta
    ((pi*exp(-ldelta + lzeta))/(sqrt(1+pi^2))),   # differentiate wrt pi
    (-sqrt(1+pi^2)*exp(-ldelta + lzeta))          # differentiate wrt ldelta
  ),
  ncol=3)

#-----------------------------------------------------------------------------
# The row vector D of the partial derivatives for beta
#-----------------------------------------------------------------------------

D.beta <- matrix(
  c(
    (pi*exp(-ldelta + lzeta)),    # differentiate wrt lzeta
    exp(-ldelta + lzeta),         # differentiate wrt pi
    (-pi*exp(-ldelta + lzeta))    # differentiate wrt ldelta
  ),
  ncol=3)

#-----------------------------------------------------------------------------
# Calculate the approximations of the variances for alpha and beta
# "sigma" denotes the 3x3 covariance matrix with log(delta) and log(zeta)
#-----------------------------------------------------------------------------

var.alpha <- D.alpha %*% sigma %*% t(D.alpha) 
var.beta <- D.beta %*% sigma %*% t(D.beta)

#-----------------------------------------------------------------------------
# The standard errors are the square roots of the variances
#-----------------------------------------------------------------------------

se.alpha <- sqrt(var.alpha)
se.beta <- sqrt(var.beta)
COOLSerdash
fonte
1
Cov(π,ζ)πζβζ/δ×π/δ=(ζπ)/δ2
1
varcov <- solve(hyperbfitalv$hessian)π,ζ,δ
3
Muito obrigado pela sua resposta, mas este é EXATAMENTE o problema, porque a parametrização deste hessian é para log (delta) e log (zeta) e não para delta e zeta! Ver o meu follow-up mensagens: stats.stackexchange.com/questions/67595/... e, especialmente, a resposta de CT Zhu aqui stats.stackexchange.com/questions/67602/...
Jen Bohold
2
é preciso obter o hessiano de pi, log (zeta), log (delta) e mu para o hessian de pi, zeta, delta e mu. Você sabe como isso pode ser feito?
Jen Bohold
2
Também tentei fazê-lo com o hessian "errado"; portanto, com log (delta) e log (zeta), depois disso, selecionei a sub-matriz e fiz os cálculos. Os resultados não estavam corretos, porque os valores eram muito grandes, mais ou menos 60.000.
Jen Bohold
-2

Possível duplicata: Erros padrão do hyperbFit?

Eu poderia apostar que algumas contas pertencem à mesma pessoa ...

YouMustWatchTheWire
fonte
Eu já vinculei esse tópico no meu post! E se você realmente ler minha pergunta, verá que NÃO QUERO usar o algoritmo de autoinicialização, mas estou perguntando sobre o método delta.
Jen Bohold