Standardfehler von hyperbolischen Verteilungsschätzungen nach der Delta-Methode?

8

Ich möchte die Standardfehler einer angepassten hyperbolischen Verteilung berechnen.

In meiner Notation ist die Dichte gegeben durch Ich verwende das HyperbolicDistr-Paket in R. Ich schätze die Parameter über den folgenden Befehl:

H.(l;;α,β,μ,δ)=α2- -β22αδK.1(δα2- -β2)exp(- -αδ2+(l- -μ)2+β(l- -μ))
hyperbFit(mydata,hessian=TRUE)

Dies gibt mir eine falsche Parametrierung. Ich ändere es mit dem hyperbChangePars(from=1,to=2,c(mu,delta,pi,zeta))Befehl in meine gewünschte Parametrierung . Dann möchte ich die Standardfehler meiner Schätzungen haben, ich kann sie für die falsche Parametrierung mit dem summaryBefehl erhalten. Dies gibt mir aber die Standardfehler für die andere Parametrierung. Laut diesem Thread muss ich die Delta-Methode verwenden (ich möchte keinen Bootstrap oder Kreuzvalidierung oder so verwenden).

Der hyperbFitCode ist hier . Und das hyperbChangeParsist hier . Daher weiß ich, dass und bleiben. Daher sind auch die Standardfehler gleich, oder?μδ

Um und in und umzuwandeln, brauche ich die Beziehung zwischen ihnen. Gemäß dem Code geschieht dies wie folgt:πζαβ

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

Wie muss ich die Delta-Methode codieren, um die gewünschten Standardfehler zu erhalten?

EDIT: Ich benutze diese Daten . Ich führe zuerst die Delta-Methode gemäß diesem Thread durch.

# fit the distribution

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

summary(hyperbfitdb) gibt folgende Ausgabe:

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 

und hyperbChangePars(from=1,to=2,hyperbfitdb$Theta)gibt die folgende Ausgabe:

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

Jetzt definiere ich die Variablen folgendermaßen:

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

Ich führe jetzt den Code aus (zweite Bearbeitung) und erhalte das folgende Ergebnis:

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

Ist das richtig? Ich wundere mich über Folgendes: Wenn ich einen Bootstrap-Algorithmus folgendermaßen verwende:

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)

Ich bekomme 119.6für sd(alpha)und 35.85für sd(beta). Die Ergebnisse sind sehr unterschiedlich? Gibt es einen Fehler oder was ist das Problem hier?

Jen Bohold
quelle

Antworten:

10

In der folgenden Lösung gehe ich hyperbPivon . Außerdem sind die in den folgenden Näherungen verwendeten Varianzen einfach die quadratischen Standardfehler, die durch after berechnet werden , also . Um die Approximation mit der Delta-Methode zu berechnen , benötigen wir die partiellen Ableitungen der Transformationsfunktion s und . Die Transformationsfunktionen für und sind gegeben durch: V a r ( X ) = S E ( X ) 2 g α ( ζ ,πsummaryhyperbFitV.einr(X.)=S.E.(X.)2g β ( ζ , π , δ ) α β g α ( ζ , π , δ )Gα(ζ,π,δ)Gβ(ζ,π,δ)αβ α

Gα(ζ,π,δ)=ζ1+π2δGβ(ζ,π,δ)=ζπδ
Die partiellen Ableitungen der Transformationsfunktion für lauten dann: Die partiellen Ableitungen der Transformationsfunktion für sind: α β
ζGα(ζ,π,δ)=1+π2δπGα(ζ,π,δ)=πζ1+π2δδGα(ζ,π,δ)=- -1+π2ζδ2
β
ζGβ(ζ,π,δ)=πδπGβ(ζ,π,δ)=ζδδGβ(ζ,π,δ)=- -πζδ2

Wenn wir die Delta-Methode auf die Transformationen anwenden, erhalten wir die folgende Näherung für die Varianz von (nehmen Sie Quadratwurzeln, um die Standardfehler zu erhalten): Die ungefähre Varianz von ist:α

V.einr(α)1+π2δ2V.einr(ζ)+π2ζ2(1+π2)δ2V.einr(π)+(1+π2)ζ2δ4V.einr(δ)+2×[πζδ2C.Öv(π,ζ)- -(1+π2)ζδ3C.Öv(δ,ζ)- -πζ2δ3C.Öv(δ,π)]]
β

V.einr(β)π2δ2V.einr(ζ)+ζ2δ2V.einr(π)+π2ζ2δ4V.einr(δ)+2×[πζδ2C.Öv(π,ζ)- -π2ζδ3C.Öv(δ,ζ)- -πζ2δ3C.Öv(π,δ)]]

Codierung in R

Der schnellste Weg, um die obigen Näherungen zu berechnen, ist die Verwendung von Matrizen. Bezeichne den Zeilenvektor, der die partiellen Ableitungen der Transformationsfunktion für oder in Bezug auf . Ferner bezeichne die Varianz-Kovarianz-Matrix von . Die Kovarianzmatrix kann durch Eingabe abgerufen werden , wo die angepasste Funktion ist. Die obige Annäherung der Varianz von ist dann Dasselbe gilt für die Approximation der Varianz vonα β ζ , π , δ Σ 3 × 3 ζ , π , δ αD.αβζ,π,δΣ3×3ζ,π,δvcov(my.hyperbFit)my.hyperbFitα

V.einr(α)D.αΣD.α
β.

In Rkann dies einfach wie folgt codiert werden:

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

Verwenden von undLog(ζ)Log(δ)

Wenn die Standardfehler / -abweichungen nur für und anstelle von und verfügbar sind , ändern sich die Transformationsfunktionen in : Die partiellen Ableitungen der Transformationsfunktion für sind dann: δ *ζ=Log(ζ)δ=Log(δ)ζδ

Gα(ζ,π,δ)=exp(ζ)1+π2exp(ζ)Gβ(ζ,π,δ)=exp(ζ)πexp(δ)
α
ζGα(ζ,π,δ)=1+π2exp(- -δ+ζ)πGα(ζ,π,δ)=πexp(- -δ+ζ)1+π2δGα(ζ,π,δ)=- -1+π2exp(- -δ+ζ)
Die partiellen Ableitungen der Transformationsfunktion für sind: β
ζGβ(ζ,π,δ)=πexp(- -δ+ζ)πGβ(ζ,π,δ)=exp(- -δ+ζ)δGβ(ζ,π,δ)=- -πexp(- -δ+ζ)
Anwenden der Delta-Methode auf Bei den Transformationen erhalten wir die folgende Näherung für die Varianz von : α
V.einr(α)(1+π2)exp(- -2δ+2ζ)V.einr(ζ)+π2exp(- -2δ+2ζ)1+π2V.einr(π)+(1+π2)exp(- -2δ+2ζ)V.einr(δ)+2×[πexp(- -2δ+2ζ)C.Öv(π,ζ)- -(1+π2)exp(- -2δ+2ζ)C.Öv(δ,ζ)- -πexp(- -2δ+2ζ)C.Öv(δ,π)]]
Die ungefähre Varianz von ist: β
V.einr(β)π2exp(- -2δ+2ζ)V.einr(ζ)+exp(- -2δ+2ζ)V.einr(π)+π2exp(- -2δ+2ζ)V.einr(δ)+2×[πexp(- -2δ+2ζ)C.Öv(π,ζ)- -π2exp(- -2δ+2ζ)C.Öv(δ,ζ)- -πexp(- -2δ+2ζ)C.Öv(δ,π)]]

Codierung in R2

Dieses Mal sigmabezeichnet die Kovarianzmatrix, jedoch einschließlich der Varianzen und Kovarianzen für und anstelle von und .ζ=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
quelle
1
@BenBohold Die Terme in den Klammern vor den Kovarianztermen sind die Produkte der jeweiligen partiellen Ableitungen der Transformationsfunktionen. Zum Beispiel: Der Term vor ist das Produkt der partiellen Ableitung wrt multipliziert mit der partiellen Ableitung wrt . Im Fall von wäre dies: . C.Öv(π,ζ)πζβζ/.δ×π/.δ=(ζπ)/.δ2
COOLSerdash
1
@ BenBohold Seltsam, aber kein Problem. Versuchen Sie, die Kovarianzmatrix folgendermaßen zu berechnen : varcov <- solve(hyperbfitalv$hessian). Funktioniert das? Danach müssen Sie die Submatrix auswählen, die nur . Am einfachsten könnte ich Ihnen helfen, wenn Sie ein voll funktionsfähiges Beispiel mit Daten bereitstellen würden (Sie müssen nicht alle Ihre Daten bereitstellen). π,ζ,δ
COOLSerdash
3
Ein großes Dankeschön für Ihre Antwort, aber das ist genau das Problem, denn die Parametrisierung dieses Hessischen ist für log (Delta) und log (Zeta) und nicht für Delta und Zeta! Siehe meine Follow-up-Beiträge: stats.stackexchange.com/questions/67595/… und insbesondere die Antwort von CT Zhu hier stats.stackexchange.com/questions/67602/…
Jen Bohold
2
man muss das Hessische von pi, log (zeta), log (delta) und mu in das Hessische von pi, zeta, delta und mu bekommen. Wissen Sie, wie das gemacht werden könnte?
Jen Bohold
2
Ich habe auch versucht, es mit dem "falschen" Hessischen zu machen, also mit log (Delta) und log (Zeta), danach habe ich die Submatrix ausgewählt und die Berechnungen durchgeführt. Die Ergebnisse waren nicht korrekt, da die Werte viel zu groß waren, also etwa 60 000.
Jen Bohold
-2

Mögliches Duplikat: Standardfehler von hyperbFit?

Ich könnte wetten, dass einige Konten derselben Person gehören ...

YouMustWatchTheWire
quelle
Ich habe diesen Thread bereits in meinem Beitrag verlinkt! Und wenn Sie meine Frage wirklich lesen, werden Sie sehen, dass ich den Bootstrap-Algorithmus NICHT verwenden möchte, aber ich frage nach der Delta-Methode.
Jen Bohold