Erreurs standard des estimations de distribution hyperbolique utilisant la méthode delta?

8

Je veux calculer les erreurs standard d'une distribution hyperbolique ajustée.

Dans ma notation, la densité est donnée par J'utilise le package HyperbolicDistr dans R. J'évalue les paramètres via la commande suivante:

H(l;α,β,μ,δ)=α2β22αδK1(δα2β2)exp(αδ2+(lμ)2+β(lμ))
hyperbFit(mydata,hessian=TRUE)

Cela me donne un mauvais paramétrage. Je le change en mon paramétrage souhaité avec la hyperbChangePars(from=1,to=2,c(mu,delta,pi,zeta))commande. Ensuite, je veux avoir les erreurs standard de mes estimations, je peux les obtenir pour le mauvais paramétrage avec la summarycommande. Mais cela me donne les erreurs standard pour l'autre paramétrage. Selon ce fil , je dois utiliser le delta méthode (je ne pas veux utiliser bootstrap ou validation croisée ou si).

Le hyperbFitcode est ici . Et le hyperbChangeParsest ici . Par conséquent, je sais que et restent les mêmes. Par conséquent, les erreurs standard sont également les mêmes, non?μδ

Pour transformer et en et j'ai besoin de la relation entre eux. Selon le code, cela se fait comme suit:πζαβ

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

Alors, comment dois-je coder la méthode delta pour obtenir les erreurs standard souhaitées?

EDIT: J'utilise ces données . J'exécute d'abord la méthode delta selon ce fil.

# fit the distribution

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

summary(hyperbfitdb) donne la sortie suivante:

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 

et hyperbChangePars(from=1,to=2,hyperbfitdb$Theta)donne la sortie suivante:

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

maintenant je définis les variables de la manière suivante:

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

Je lance maintenant le code (deuxième édition) et j'obtiens le résultat suivant:

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

Est-ce correct? Je me pose des questions sur ce qui suit: si j'utilise un algorithme d'amorçage de la manière suivante:

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)

Je reçois 119.6pour sd(alpha)et 35.85pour sd(beta). Les résultats sont très différents? Y a-t-il une erreur ou quel est le problème ici?

Jen Bohold
la source

Réponses:

10

Dans la solution suivante, je suppose hyperbPiêtre . De plus, les variances utilisées dans les approximations ci-dessous sont simplement les erreurs standard au carré calculées par après , donc . Pour calculer l'approximation à l'aide de la méthode delta , nous avons besoin des dérivées partielles de la fonction de transformation s et . Les fonctions de transformation pour et sont données par: πsummaryhyperbFitVar(X)=SE(X)2gα(ζ,π,δ)gβ(ζ,π,δ)αβ

gα(ζ,π,δ)=ζ1+π2δgβ(ζ,π,δ)=ζπδ
Les dérivées partielles de la fonction de transformation pour sont alors: Les dérivées partielles de la fonction de transformation pour sont: α
ζgα(ζ,π,δ)=1+π2δπgα(ζ,π,δ)=πζ1+π2δδgα(ζ,π,δ)=1+π2ζδ2
β
ζgβ(ζ,π,δ)=πδπgβ(ζ,π,δ)=ζδδgβ(ζ,π,δ)=πζδ2

En appliquant la méthode delta aux transformations, nous obtenons l'approximation suivante pour la variance de (prendre des racines carrées pour obtenir les erreurs standard): La variance approximative de est:α

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(π,δ)]

Codage dans R

Le moyen le plus rapide pour calculer les approximations ci-dessus est d'utiliser des matrices. Notons le vecteur ligne contenant les dérivées partielles de la fonction de transformation pour ou par rapport à . De plus, notons la matrice variance-covariance de . La matrice de covariance peut être récupérée en tapant où est la fonction ajustée. L'approximation ci-dessus de la variance de est alors Il en va de même pour l'approximation de la variance deDαβζ,π,δΣ3×3ζ,π,δvcov(my.hyperbFit)my.hyperbFitα

Var(α)DαΣDα
β.

Dans R, cela peut être facilement codé comme ceci:

#-----------------------------------------------------------------------------
# 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)

Utilisation de etlog(ζ)log(δ)

Si les erreurs / écarts standard ne sont disponibles que pour et au lieu de et , les fonctions de transformation passent à : Les dérivées partielles de la fonction de transformation pour sont alors: ζ=log(ζ)δ=log(δ)ζδ

gα(ζ,π,δ)=exp(ζ)1+π2exp(ζ)gβ(ζ,π,δ)=exp(ζ)πexp(δ)
α
ζgα(ζ,π,δ)=1+π2exp(δ+ζ)πgα(ζ,π,δ)=πexp(δ+ζ)1+π2δgα(ζ,π,δ)=1+π2exp(δ+ζ)
Les dérivées partielles de la fonction de transformation pour sont: β
ζgβ(ζ,π,δ)=πexp(δ+ζ)πgβ(ζ,π,δ)=exp(δ+ζ)δgβ(ζ,π,δ)=πexp(δ+ζ)
Application de la méthode delta à les transformations, on obtient l'approximation suivante pour la variance de : α
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(δ,π)]
La variance approximative de est: β
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(δ,π)]

Codage en R2

Cette fois, sigmadénote la matrice de covariance mais incluant les variances et covariances pour et au lieu de et .ζ=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
la source
1
@BenBohold Les termes entre parenthèses avant les termes de covariance sont les produits des dérivées partielles respectives des fonctions de transformation. Par exemple: Le terme avant est le produit de la dérivée partielle wrt multipliée par la dérivée partielle wrt . Dans le cas de ce serait: . Cov(π,ζ)πζβζ/δ×π/δ=(ζπ)/δ2
COOLSerdash
1
@BenBohold Strange, mais pas de problème. Essayez de calculer la matrice de covariance de cette façon: varcov <- solve(hyperbfitalv$hessian). Est-ce que ça marche? Après cela, vous devrez sélectionner la sous-matrice contenant uniquement . La façon la plus simple dont je pourrais vous aider serait de fournir un exemple pleinement fonctionnel avec des données (vous n'avez pas à fournir toutes vos données). π,ζ,δ
COOLSerdash
3
Un grand merci pour votre réponse, mais c'est EXACTEMENT le problème, car le paramétrage de cette toile de jute est pour log (delta) et log (zeta) et non pour delta et zeta! Voir mes articles de suivi: stats.stackexchange.com/questions/67595/… et surtout la réponse de CT Zhu ici stats.stackexchange.com/questions/67602/…
Jen Bohold
2
il faut obtenir la toile de jute de pi, log (zêta), log (delta) et mu dans la toile de jute de pi, zeta, delta et mu. Savez-vous comment cela pourrait se faire?
Jen Bohold
2
J'ai également essayé de le faire avec la "mauvaise" toile de jute, donc avec log (delta) et log (zeta), après cela, j'ai sélectionné la sous-matrice et fait les calculs. Les résultats n'étaient pas corrects, car les valeurs étaient beaucoup trop grandes, comme environ 60 000 environ.
Jen Bohold
-2

Duplicata possible: erreurs standard d'hyperbFit?

Je pourrais parier que certains comptes appartiennent à la même personne ...

YouMustWatchTheWire
la source
J'ai déjà lié ce fil dans mon post! Et si vous lisez vraiment ma question, vous verrez que je NE VEUX PAS utiliser l'algorithme de bootstrap, mais je pose des questions sur la méthode delta.
Jen Bohold