*! extreme 1.2.0 17 May 2017
*! Copyright (C) 2015-17 David Roodman
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
* You should have received a copy of the GNU General Public License
* along with this program. If not, see .
mata
mata clear
mata set matastrict on
mata set mataoptimize on
mata set matalnum off
// similar to selectindex(), introduced in Stata 13
real vector extreme_selectindex(real vector v) {
real scalar rows; real vector retval
rows = rows(v)
retval = select(rows>1? 1::rows : 1..cols(v), v)
return(cols(retval)? retval : J(0,1,0))
}
real matrix extreme_clone(real matrix X) {
real matrix Y
return(Y = X)
}
// Mata interface for estimator: one-time set-up. Returns moptimize() estimation object
transmorphic extreme_init(real scalar N, real scalar GEV, | string scalar wttype, real colvector wt) {
transmorphic M
M = moptimize_init()
moptimize_init_evaluator(M, &extreme_lf2())
moptimize_init_evaluatortype(M, "lf2")
moptimize_init_depvar(M, 1, J(N,0,0)) // cue moptimize() about sample size
if (GEV) moptimize_init_eq_name(M, 1, "mu")
moptimize_init_eq_name(M, 1+GEV, "lnsig")
moptimize_init_eq_name(M, 2+GEV, "xi")
moptimize_init_search(M, "off")
moptimize_init_userinfo(M, 1, GEV)
if (wttype!="") {
moptimize_init_weight(M, wt)
moptimize_init_weighttype(M, wttype)
}
return (M)
}
// shift available order stats so lowest in rightmost col, for programming convenience in extreme_lf2()
// modifies argument and returns the max order available by observation (scalar if constant)
real matrix extreme_est_prep(real matrix Z) {
real colvector r; real scalar i
if (cols(Z)>1 & hasmissing(Z)) {
r = cols(Z) :- rowmax((Z:==.) :* (cols(Z)..1))
if (isview(Z)) Z = extreme_clone(Z)
for (i=rows(Z); i; i--)
if (r[i]1? extreme_selectindex(abs(xi):<1e-10) : (abs(xi):<1e-10? 1::rows(x) : J(0,0,0) )
if (rows(gumbel)) xgumbel=x[gumbel,]
inv_xi = 1:/xi
lny = ln(y = 1:+xi:*x)
if (GEV) {
f = (R>1? y[,R] : y):^-inv_xi
if (rows(gumbel)) f[gumbel] = exp(-(R>1? xgumbel[,R] : xgumbel))
}
if (todo) {
real colvector neg_lny_x, neg_lny_xi, lnf_x, lnf_xi, lng_x, lng_xi, lng_lnsig, inv_sig, neg_inv_sig
lnf_x = (-1):/y
neg_lny_x = xi:*lnf_x
neg_lny_xi = x :*lnf_x
lnf_xi = inv_xi:*(inv_xi:*lny+neg_lny_xi)
if (rows(gumbel)) lnf_xi[gumbel,] = 0.5*xgumbel:*xgumbel
lng_x = lnf_x + neg_lny_x
lng_xi = lnf_xi + neg_lny_xi; if (R>1) lng_xi = quadrowsum(lng_xi)
lng_lnsig = x :* lng_x ; if (R>1) lng_lnsig = quadrowsum(lng_lnsig)
if (GEV) {
f_x = f :* (R>1? lnf_x [,R] : lnf_x )
lng_xi = lng_xi - f :* (R>1? lnf_xi[,R] : lnf_xi)
neg_inv_sig = -(inv_sig = 1:/sig)
lng_mu = neg_inv_sig :* ((R>1? quadrowsum(lng_x):lng_x) - f_x)
lng_lnsig = lng_lnsig - (R>1? x[,R] : x) :*f_x
}
S = -r:-lng_lnsig, lng_xi; if (GEV) S = lng_mu, S
if (todo==2) {
real colvector neg_lny_xx, neg_lny_xixi, lnf_xx, lnf_xxi, lnf_xixi, lng_xx, lng_xxi, lng_xixi, sig_lng_mulnsig, lny_xxi, t
real matrix d2lnsiglnsig, d2lnsigxi, d2xixi, d2mumu, d2mulnsig, d2muxi
real scalar v; v = 1
neg_lny_xx = neg_lny_x:*neg_lny_x
lny_xxi = lnf_x:*lnf_x
neg_lny_xixi = neg_lny_xi:*neg_lny_xi
lnf_xx = xi:*lny_xxi
lnf_xxi = x :*lny_xxi
lnf_xixi = inv_xi:*(neg_lny_xixi-2*lnf_xi)
if (rows(gumbel)) lnf_xixi[gumbel,] = (-1.3333333333333333)*lnf_xi[gumbel,]:*xgumbel
lng_xx = lnf_xx + neg_lny_xx
lng_xxi = lnf_xxi - lny_xxi
lng_xixi = lnf_xixi + neg_lny_xixi
if (GEV)
if (R==1) {
lng_xx = lng_xx - f:*(lnf_xx + lny_xxi ) // lny_xxi = lnf_x^2
lng_xxi = lng_xxi - f:*(lnf_xxi + lnf_xi :*lnf_x )
lng_xixi = lng_xixi - f:*(lnf_xixi + lnf_xi :*lnf_xi )
} else {
lng_xx [,R] = lng_xx [,R] - f:*(lnf_xx [,R] + lny_xxi[,R] )
lng_xxi [,R] = lng_xxi [,R] - f:*(lnf_xxi [,R] + (t=lnf_xi [,R]):*lnf_x[,R])
lng_xixi[,R] = lng_xixi[,R] - f:*(lnf_xixi[,R] + t:*t )
}
sig_lng_mulnsig = lng_x + x:*lng_xx
if (GEV)
if (R==1)
sig_lng_mulnsig = sig_lng_mulnsig - f_x
else
sig_lng_mulnsig[,R] = sig_lng_mulnsig[,R] - f_x
d2lnsigxi = -moptimize_util_matsum(M, 1+GEV, 2+GEV, quadrowsum(x :* lng_xxi ), v)
d2lnsiglnsig = moptimize_util_matsum(M, 1+GEV, 1+GEV, quadrowsum(x :* sig_lng_mulnsig), v)
d2xixi = moptimize_util_matsum(M, 2+GEV, 2+GEV, quadrowsum(lng_xixi ), v)
if (GEV) {
if (rows(sig)==1) {
d2mumu = moptimize_util_matsum(M, 1, 1, quadrowsum(lng_xx ), v)/(sig*sig)
d2mulnsig = moptimize_util_matsum(M, 1, 2, quadrowsum(sig_lng_mulnsig), v)/sig
d2muxi = moptimize_util_matsum(M, 1, 3, quadrowsum(lng_xxi ), v)/(-sig)
} else {
d2mumu = moptimize_util_matsum(M, 1, 1, quadrowsum(lng_xx :*neg_inv_sig:*neg_inv_sig), v)
d2mulnsig = moptimize_util_matsum(M, 1, 2, quadrowsum(sig_lng_mulnsig :*inv_sig ), v)
d2muxi = moptimize_util_matsum(M, 1, 3, quadrowsum(lng_xxi :*neg_inv_sig ), v)
}
H = d2mumu,d2mulnsig,d2muxi \ d2mulnsig',d2lnsiglnsig,d2lnsigxi \ d2muxi',d2lnsigxi',d2xixi
} else
H = d2lnsiglnsig,d2lnsigxi \ d2lnsigxi',d2xixi
}
}
if (R>1) {
_editmissing(x, 0) // should only be missing in cells representing missing statistics of order>1. Editing to 0 results in 0 contribution to likelihood.
_editmissing(xgumbel, 0) // ...or if sig overflowed, in which case lng will become missing below
lny = ln(1:+xi:*x)
}
lng = quadrowsum((inv_xi:+1):*lny, 1)
if (rows(gumbel)) lng[gumbel] = quadrowsum(xgumbel)
lng = (R>1? -r*ln(sig) : -ln(sig)) :- lng
if (GEV)
lng = lng - f
}
// parameterize GEV or GPD by (mu,) z_p, xi, where z_p=return level for probability p
function extreme_return_lf(transmorphic M, real rowvector b, real colvector lng) {
real matrix x, y; real scalar zp, mu, sig, xi, GEV, r, yp, gumbel
GEV = moptimize_util_userinfo(M, 1)
xi = moptimize_util_xb(M, b, 2+GEV)
r = moptimize_util_userinfo(M, 3)
zp = moptimize_util_xb(M, b, 1+GEV)
yp = moptimize_util_userinfo(M, 4)
mu = GEV? moptimize_util_xb(M, b, 1) : 0
gumbel = abs(xi)<1e-10
sig = (zp - mu) / (gumbel? ln(yp) : (yp^xi-1)/xi)
x = moptimize_util_userinfo(M, 2); x = (GEV? (x :- mu) : x) / sig
lng = -(cols(x)>1? r*ln(sig) : ln(sig)) :- (gumbel? x : (1+1/xi)*ln(y = 1:+xi*x))
if (GEV) lng = lng - (gumbel? exp(-x) : y:^-(1/xi))
}
// parameterize GEV or GPD by (mu,) p, xi, p=return rate
function extreme_rate_lf(transmorphic M, real rowvector b, real colvector lng) {
real matrix x, y; real scalar p, mu, sig, xi, GEV, r, yp, zp, gumbel
GEV = moptimize_util_userinfo(M, 1)
xi = moptimize_util_xb(M, b, 2+GEV)
r = moptimize_util_userinfo(M, 3)
p = moptimize_util_xb(M, b, 1+GEV)
zp = moptimize_util_userinfo(M, 4) // extremal magnitude of interest
if (GEV) {
mu = moptimize_util_xb(M, b, 1)
yp = (-1):/ln(1:-p)
} else {
mu = 0
yp = 1:/p
}
gumbel = abs(xi)<1e-10
sig = (zp :- mu) :/ (gumbel? ln(yp) : (yp:^xi:-1):/xi)
x = moptimize_util_userinfo(M, 2); x = (GEV? (x :- mu) : x) :/ sig
lng = -(cols(x)>1? r*ln(sig) : ln(sig)) :- (gumbel? x : (1:+1:/xi):*ln(y = 1:+xi:*x))
if (GEV) lng = lng - (gumbel? exp(-x) : y:^-(1:/xi))
}
// parameterize GEV or GPD by (mu,) z_p, xi, where z_p=return level for probability p
function extreme_return_lf1(transmorphic M, real scalar todo, real rowvector b, real colvector lng, real matrix S, real matrix H) {
real matrix x, lny; real colvector lng_mu, f, f_x, r
real scalar GEV, R, yp, lnyp, zp, sig, xi, mu, zpmmu, y, gumbel, inv_xi
pragma unused H
GEV = moptimize_util_userinfo(M, 1)
zp = moptimize_util_xb(M, b, 1+GEV)
xi = moptimize_util_xb(M, b, 2+GEV)
mu = GEV? moptimize_util_xb(M, b, 1) : 0
lnyp = ln(yp = moptimize_util_userinfo(M, 4))
gumbel = abs(xi):<1e-10
sig = gumbel? lnyp : (yp^xi-1) / xi; sig = (zpmmu = zp - mu) / sig
if (GEV) {
x = moptimize_util_userinfo(M, 2) :- mu
R = cols(x)
x = x / sig
r = moptimize_util_userinfo(M, 3)
} else {
x = moptimize_util_userinfo(M, 2) / sig
r = R = 1
}
inv_xi = 1/xi
lny = ln(y = 1:+xi*x)
if (GEV)
f = gumbel? exp(-(R>1? x[,R] : x)) : (R>1? y[,R] : y):^-inv_xi
if (todo) {
real colvector neg_lny_x, neg_lny_xi, lnf_x, lnf_xi, lng_x, lng_xi, lng_lnsig, neg_inv_sig
lnf_x = (-1):/y
neg_lny_x = xi *lnf_x
neg_lny_xi = x :*lnf_x
lnf_xi = gumbel? 0.5*x:*x : inv_xi*(inv_xi*lny+neg_lny_xi)
lng_x = lnf_x + neg_lny_x
lng_xi = lnf_xi + neg_lny_xi; if (R>1) lng_xi = rowsum(lng_xi)
lng_lnsig = x :* lng_x ; if (R>1) lng_lnsig = rowsum(lng_lnsig)
if (GEV) {
f_x = f :* (R>1? lnf_x [,R] : lnf_x )
lng_xi = lng_xi - f :* (R>1? lnf_xi[,R] : lnf_xi)
neg_inv_sig = (-1):/sig
lng_mu = neg_inv_sig :* ((R>1? rowsum(lng_x):lng_x) - f_x)
lng_lnsig = lng_lnsig - (R>1? x[,R] : x) :*f_x
}
lng_lnsig = -r :- lng_lnsig
S = lng_lnsig / zpmmu, lng_lnsig * (gumbel? -.5*lnyp*lnyp : 1/xi-lnyp/(1-yp^(-xi))) + lng_xi; if (GEV) S = lng_mu - S[,1] , S // change of variables from mu,lnsig,xi to mu,zp,xi
}
if (R>1) {
_editmissing(x, 0) // should only be missing in cells representing missing statistics of order>1. Editing to 0 results in 0 contribution to likelihood.
if (gumbel) _editmissing(x, 0) // ...or if sig overflowed, in which case lng will become missing below
lny = ln(1:+xi*x)
}
lng = gumbel? rowsum(x) : (inv_xi+1) * rowsum(lny, 1)
lng = -(R>1? r*ln(sig) : ln(sig)) :- lng
if (GEV)
lng = lng - f
}
// parameterize GEV or GPD by (mu,) p, xi, p=return rate
function extreme_rate_lf1(transmorphic M, real scalar todo, real rowvector b, real colvector lng, real matrix S, real matrix H) {
real matrix x, lny; real colvector f, f_x, r
real scalar GEV, R, yp, lnyp, zp, sig, xi, mu, zpmmu, y, gumbel, inv_xi, p, fp
pragma unused H
GEV = moptimize_util_userinfo(M, 1)
xi = moptimize_util_xb(M, b, 2+GEV)
p = moptimize_util_xb(M, b, 1+GEV)
zp = moptimize_util_userinfo(M, 4) // extremal magnitude of interest
if (GEV) {
mu = moptimize_util_xb(M, b, 1)
yp = -1/ln(1-p)
} else {
mu = 0
yp = 1/p
}
lnyp = ln(yp)
gumbel = abs(xi):<1e-10
sig = gumbel? lnyp : (yp^xi-1) / xi; sig = (zpmmu = zp - mu) / sig
if (GEV) {
mu = moptimize_util_xb(M, b, 1)
x = moptimize_util_userinfo(M, 2) :- mu
R = cols(x)
x = x / sig
r = moptimize_util_userinfo(M, 3)
} else {
x = moptimize_util_userinfo(M, 2) / sig
r = R = 1
}
inv_xi = 1/xi
lny = ln(y = 1:+xi*x)
if (GEV)
f = gumbel? exp(-(R>1? x[,R] : x)) : (R>1? y[,R] : y):^-inv_xi
if (todo) {
real colvector neg_lny_x, neg_lny_xi, lnf_x, lnf_xi, lng_x, lng_xi, lng_lnsig, neg_inv_sig, lng_mu, lng_p, xp
lnf_x = (-1):/y
neg_lny_x = xi *lnf_x
neg_lny_xi = x :*lnf_x
lnf_xi = gumbel? 0.5*x:*x : inv_xi*(inv_xi*lny+neg_lny_xi)
lng_x = lnf_x + neg_lny_x
lng_xi = lnf_xi + neg_lny_xi; if (R>1) lng_xi = rowsum(lng_xi)
lng_lnsig = x :* lng_x ; if (R>1) lng_lnsig = rowsum(lng_lnsig)
if (GEV) {
f_x = f :* (R>1? lnf_x [,R] : lnf_x )
lng_xi = lng_xi - f :* (R>1? lnf_xi[,R] : lnf_xi)
neg_inv_sig = (-1):/sig
lng_mu = neg_inv_sig :* ((R>1? rowsum(lng_x):lng_x) - f_x)
lng_lnsig = lng_lnsig - (R>1? x[,R] : x) :*f_x
}
lng_lnsig = -r :- lng_lnsig
// change of variables from mu,lnsig,xi to mu,p,xi
lnyp = ln(yp = 1 + xi * (xp = zpmmu/sig)) // not conceptually the same "yp" as above--here, y at the quantile defined by p
fp = yp^(-inv_xi)
lng_p = lng_lnsig / (zpmmu/sig*fp/yp); if (GEV) lng_p = lng_p / exp(-fp)
lng_xi = lng_xi + lng_lnsig * (inv_xi - inv_xi*inv_xi / xp * yp * lnyp)
S = lng_p, lng_xi
if (GEV) S = lng_mu - lng_lnsig/zpmmu, S
}
if (R>1) {
_editmissing(x, 0) // should only be missing in cells representing missing statistics of order>1. Editing to 0 results in 0 contribution to likelihood.
if (gumbel) _editmissing(x, 0) // ...or if sig overflowed, in which case lng will become missing below
lny = ln(1:+xi*x)
}
lng = gumbel? rowsum(x) : (inv_xi+1) * rowsum(lny, 1)
lng = -(R>1? r*ln(sig) : ln(sig)) :- lng
if (GEV)
lng = lng - f
}
// represents sum of terms of form C * xi^a * sig^b * u^l * (ln u)^d * x^n / y^m
class poly {
real scalar len // number of terms
real matrix coef // one col for each poly, one row (if more than one) for each obs
real matrix powers // one row for each term's a, b, l, d, n, m 6-tuple
static pointer (real colvector function) scalar pfnPi // expectation function to be used for all polys
static pointer (real colvector) scalar pxi, psig, pwt // xi & sig values, observation weights
static real colvector gumbel // rows of xi near zero
static real scalar tol // tolerance defining "near"
static real scalar N // number of obs represented in calculations--1 if no covariates
static class poly scalar term()
void _C(), _xi(), _sig(), debug()
class poly scalar C(), x(), u(), lnu(), dxi(), sig()
static class poly scalar sum()
static class poly rowvector times()
void set_xi_sig_wt_fnPi()
void cleanup(), equals()
real matrix E()
protected real colvector _E()
void new()
}
void poly::new() {
len = 0
coef = J(0,1,0)
powers = J(0,6,0)
if (pfnPi==NULL) pfnPi = &Pi_GEV()
if (pwt==NULL) pwt = pxi = &J(0,1,0)
}
void poly::debug() {
"coef,powers"
coef,powers
"gumbel"
gumbel
}
void poly::set_xi_sig_wt_fnPi(real colvector xi, real scalar tol, real colvector sig, real colvector wt, pointer (real colvector function) scalar pPi) {
pxi = ξ pwt = &wt; psig = &sig; pfnPi = pPi; this.tol = tol
N = max((rows(xi), rows(wt), rows(sig)))
gumbel = extreme_selectindex(abs(xi):1) {
p = order(powers, -1..-6) // put high negative powers of xi at bottom so backward loop in E() evaluates most potentially explosive terms first
coef = coef[p,]
_collate(powers,p)
p = 1 \ rowsum(powers[|2,.\.,.|] :!= powers[|.,.\rows(coef)-1,.|]):!=0 // rows for new term types
powers = powers[extreme_selectindex(p),]
coef = designmatrix(runningsum(p)) ' coef
}
p = extreme_selectindex(rowsum(coef:!=0):>0)
coef = coef[p,]
powers = powers[p,]
}
len = rows(coef)
}
// differentiate a poly wrt xi--in sense that E[result]=dE[argument]/dxi
class poly scalar poly::dxi() {
class poly scalar t1, t2, t3, t4, retval
t1 = this
t1.coef = t1.coef :* -t1.powers[,5] // -n * u^l * (ln u)^d * x^n / y^m
t2 = t1.lnu().x(-1) // // -n * u^l * (ln u)^(d+1) * x^(n-1) / y^m
retval = poly::sum((t1,t2)); retval._xi(-1) // (t1+t2)/xi
t3 = this.lnu(); t3.coef = t3.coef :* (t3.powers[,6]-t3.powers[,5]) // (m-n) * [u^l * (ln u)^(d+1) * x^n / y^m]
t4 = this; t4.coef = coef :* powers[,1]; t4.powers[,1] = t4.powers[,1] :- 1 // Dxi^a = a*xi^(a-1)
retval = poly::sum((retval, t3, t4)) // dxipoly_dxi * [u^l*(ln u)^d*x^n/y^m] + xipoly * {-n/xi*([u^l*(ln u)^d*x^n/y^m] + [u^l*(ln u)^(d+1)*x^(n-1)/y^m]) + (m-n)*[u^l*(ln u)^(d+1)*x^n/y^m] }
retval.cleanup()
return(retval)
}
// expectations of a poly, optionally restricted to observations p. gumbel holds indexes of small xi *after* any restriction by p
real colvector poly::_E(real colvector xi, | real colvector p) {
real scalar i; real colvector E
pointer (real matrix) scalar _psig; pointer (real colvector) _pxi
_psig = rows(p) & rows(*psig)>1? &((*psig)[p,]) : psig
_pxi = rows(p) & rows(xi )>1? &( xi [p,]) : &xi
E = J(N, 1, 0)
for (i=len; i; i--)
E = E + quadrowsum((*_pxi):^powers[i,1] :* (*_psig):^powers[i,2] :* coef[i] :* (-1):^powers[i,5] :* BackDiff(powers[i,3],powers[i,4],powers[i,5],powers[i,6], xi, pfnPi))
return (E)
}
// weighted, observation-wise expectations of a rowvector of polys
real matrix poly::E(class poly rowvector p) {
real scalar i; real matrix retval
retval = J(N, cols(p), 0)
for (i=cols(p); i; i--) {
if (rows(*pxi)==rows(gumbel))
retval[,i] = ((tol:-*pxi):* p[i]._E(-tol, gumbel) + (tol:+*pxi):* p[i]._E(tol, gumbel))/(2*tol)
else {
retval[,i] = p[i]._E(*pxi)
if (rows(gumbel)) retval[gumbel,i] = ((tol:-(*pxi)[gumbel]):* p[i]._E(-tol, gumbel) + (tol:+(*pxi)[gumbel]):* p[i]._E(tol, gumbel))/(2*tol)
}
}
return(retval :* *pwt)
}
// hold pre-computed values of gamma(), digamma(), and trigamma() of a+b*xi, a and b certain integers, for speed
class gammastore {
static transmorphic _gamma, _digamma, _trigamma
static real scalar amin, amax, bmin_gamma, bmin_digamma, bmin_trigamma, bmax_gamma, bmax_digamma, bmax_trigamma
static pointer(real colvector) scalar pxi
void init()
real matrix gamma(), digamma(), trigamma()
}
// precompute gamma(), digamma(), trigamma() of a+b*xi, for pre-set ranges of a,b and given colvector xi, as well a+b*(+/- tol)
// associative array indexed by a and i, where i=-1 means -tol, i=0 means xi, and i=1 means +tol
void gammastore::init(real colvector xi, real scalar tol, real scalar R) {
real scalar a, i; pointer(real colvector) rowvector p
amin = 1
amax = 9+R
bmin_gamma = 0
bmin_digamma = 0
bmin_trigamma = 0
bmax_gamma = 3
bmax_digamma = 2
bmax_trigamma = 1
pxi = &xi
_gamma = asarray_create("real", 2)
_digamma = asarray_create("real", 2)
_trigamma = asarray_create("real", 2)
p = &(-tol), &xi, &tol
for (i=3; i; i--)
for (a=amin; a<=amax; a++) {
asarray( _gamma, (a,i-2), ::gamma (a:+*p[i]*( bmin_gamma.. bmax_gamma)))
asarray( _digamma, (a,i-2), ::digamma (a:+*p[i]*( bmin_digamma.. bmax_digamma)))
asarray(_trigamma, (a,i-2), ::trigamma(a:+*p[i]*(bmin_trigamma..bmax_trigamma)))
}
}
real matrix gammastore::gamma (real scalar a, real rowvector b, real colvector xi) return (asarray(_gamma, (a,(&xi==pxi? 0 : sign(xi))))[b==0?1:.,b:-(bmin_gamma -1)])
real matrix gammastore::digamma (real scalar a, real rowvector b, real colvector xi) return (asarray(_digamma, (a,(&xi==pxi? 0 : sign(xi))))[b==0?1:.,b:-(bmin_digamma -1)])
real matrix gammastore::trigamma(real scalar a, real rowvector b, real colvector xi) return (asarray(_trigamma, (a,(&xi==pxi? 0 : sign(xi))))[b==0?1:.,b:-(bmin_trigamma-1)])
// for positive integers
real scalar tetragamma(x) return (x==1? -2.4041138063191885707995 : -2.4041138063191885707995 + 2*sum((1::x-1):^-3))
//real scalar polygamma(m, x) return (-(-1)^m*factorial(m)*sum((x..x+100000):^(-m-1)))
// dth derivative of 1/(.) at a+bx
real matrix Pi_GPD(real scalar d, real scalar a, real rowvector b, real colvector x)
return ((-1)^d*factorial(d)*(1:/(a:+x*b)):^(d+1))
// nth derivative of gamma(), not lngamma(), at a+bx, where x prestored in gammastore class. For n>2, only works for a=1,2,b=0, in which case returns a scalar.
// xi=+/-1: evaluate at xi=+/-tol. xi==0: evaluate at preset xi
real matrix Pi_GEV(real scalar d, real scalar a, real rowvector b, real colvector xi) return(dgamma(d, a, b, xi))
real matrix dgamma(real scalar d, real scalar a, real rowvector b, real colvector xi) {
class gammastore scalar _gammastore
if (!d ) return(_gammastore.gamma(a,b,xi))
if (d==1) return(dgamma(0,a,b,xi):*_gammastore.digamma(a,b,xi))
if (d==2) return(dgamma(1,a,b,xi):*_gammastore.digamma(a,b,xi) + dgamma(0,a,b,xi):*_gammastore.trigamma(a,b,xi))
return(dgamma(2,a,b,xi):*_gammastore.digamma(a,b,xi) + 2*dgamma(1,a,b,xi):*_gammastore.trigamma(a,b,xi) + dgamma(0,a,b,xi)*tetragamma(a))
}
// compute 1/xi^n * backdiff_n of Pi^(d)(1+l+m*xi)
real colvector BackDiff(real scalar l, real scalar d, real scalar n, real scalar m, real colvector xi, pointer (real scalar function) pfnPi) {
real rowvector k
if (n) {
k = 0..n
return (quadrowsum((-1):^k:*comb(n,k):*(*pfnPi)(d,1+l,m:-k,xi)):/xi:^n)
}
return ((*pfnPi)(d,1+l,m,xi))
}
// do the bulk of the Cox-Snell bias-correction work, returning A & K as a pointer pair
pointer(real matrix) rowvector extreme_CSByR(real scalar GEV, real scalar R, real colvector lnsig, real colvector xi, real colvector wt, real matrix X_mu, real matrix X_lnsig, real matrix X_xi, real scalar X_mu_1, real scalar X_lnsig_1, real scalar X_xi_1) {
class poly scalar E, lny_x, lnu_x, lnu_xi, u_x, lng_x, lny_xx, lny_xxi, lny_xixi, lnu_xx, lnu_xxi, lnu_xixi, u_xx, u_xxi, u_xixi, lng_xx, lng_xxi, lng_xixi, lny_xxx, lny_xxxi, lny_xxixi,
lny_xixixi, lnu_xxx, lnu_xxxi, lnu_xxixi, lnu_xixixi, u_xxx, u_xxxi, u_xxixi, u_xixixi, lng_xxx, lng_xxxi, lng_xxixi, lng_xixixi, xlng_x, xlng_xx, x2lng_xx, xlng_xxi, lng_mumu, lng_mulnsig,
lng_lnsiglnsig, lng_muxi, lng_lnsigxi, lng_mumumu, lng_mumulnsig, lng_mulnsiglnsig, lng_lnsiglnsiglnsig, lng_mumuxi, lng_mulnsigxi,
lng_muxixi, lng_lnsiglnsigxi, lng_lnsigxixi, lng_xx_xi, lng_xxi_xi, xlng_xxi_xi, u_xxi_xi, u_xixi_xi, xu_xxi_xi,
lng_xixi_xi, lng_mumu_xi, lng_mulnsig_xi, lng_muxi_xi, lng_lnsigxi_xi
class poly scalar lng_lnsiglnsig_xi
class poly rowvector D2lng, D3lng, D2lng_xi
real rowvector E_D2lng, E_D3lng
real scalar sig, i, k
real matrix K, A, _E_D2lng, _E_D3lng, D, L, Q, QQ
pointer (real scalar function) pfnPi
class gammastore scalar _gammastore
sig = exp(lnsig)
pfnPi = GEV? &Pi_GEV() : &Pi_GPD()
if (rows(xi)==1 & rows(lnsig)==1) wt = sum(wt)
E.set_xi_sig_wt_fnPi(xi, .01, sig, wt, pfnPi)
if (GEV) _gammastore.init(xi, .01, R)
lny_x = term( 1,1,0,0,0,1) // xi/y
lnu_x = term(-1,0,0,0,0,1) // -1/y
lnu_xi = poly::sum((term(-1,-1,0,1,0,0),term(-1,-1,0,0,1,1))) // -1/xi*x/y-1/xi*ln(u)
lng_x = poly::sum((lny_x. C(-1),lnu_x)) // -lny_x + lnu_x
lny_xx = term(-1,2,0,0,0,2) // -xi^2/y^2
lny_xxi = term(1,0,0,0,0,2) // 1/y^2
lny_xixi = term(-1,0,0,0,2,2) //-x^2/y^2
lnu_xx = term(1,1,0,0,0,2) // xi/y^2
lnu_xxi = term(1,0,0,0,1,2) // x/y^2
lnu_xixi = poly::sum((term(1,-1,0,0,2,2),term(2,-2,0,1,0,0),term(2,-2,0,0,1,1))) // 1/xi*x^2/y^2+2/xi^2*ln(u)+2/xi^2*x/y
lng_xx = poly::sum((lny_xx. C(-1),lnu_xx )) // -lny_xx + lnu_xx
lng_xxi = poly::sum((lny_xxi. C(-1),lnu_xxi )) // -lny_xxi + lnu_xxi
lng_xixi = poly::sum((lny_xixi.C(-1),lnu_xixi)) // -lny_xixi + lnu_xixi
lny_xxx = term( 2,3,0,0,0,3) // 2*xi^3/y^3
lny_xxxi = term(-2,1,0,0,0,3) // -2*xi/y^3
lny_xxixi = term(-2,0,0,0,1,3) // -2*x/y^3
lny_xixixi = term( 2,0,0,0,3,3) // 2*x^3/y^3
lnu_xxx = term(-2,2,0,0,0,3) // -2*xi^2/y^3
lnu_xxxi = poly::sum((term(2,0,0,0,0,3),term(-1,0,0,0,0,2))) // 2/y^3-1/y^2
lnu_xxixi = term(-2,0,0,0,2,3) // -2*x^2/y^3
lnu_xixixi = poly::sum((term(-6,-3,0,0,1,1),term(-3,-2,0,0,2,2),term(-2,-1,0,0,3,3),term(-6,-3,0,1,0,0))) // -6/xi^3*x/y-3/xi^2*x^2/y^2-2/xi*x^3/y^3-6/xi^3*ln(u)
lng_xxx = poly::sum((lny_xxx. C(-1),lnu_xxx )) // -lny_xxx +lnu_xxx
lng_xxxi = poly::sum((lny_xxxi. C(-1),lnu_xxxi )) // -lny_xxxi +lnu_xxxi
lng_xxixi = poly::sum((lny_xxixi. C(-1),lnu_xxixi )) // -lny_xxixi +lnu_xxixi
lng_xixixi = poly::sum((lny_xixixi.C(-1),lnu_xixixi)) // -lny_xixixi +lnu_xixixi
if (GEV) {
class poly scalar lnu_x_lnu_x, lnu_x_lnu_xi, lnu_xi_lnu_xi
lnu_x_lnu_x = poly::times(lnu_x, lnu_x)
lnu_x_lnu_xi = poly::times(lnu_x, lnu_xi)
lnu_xi_lnu_xi = poly::times(lnu_xi, lnu_xi)
u_x = lnu_x.u() // u:*lnu_x
u_xx = (poly::sum((lnu_x_lnu_x , lnu_xx ))).u() // u:*lnu_x :*lnu_x + u:*lnu_xx
u_xxi = (poly::sum((lnu_x_lnu_xi , lnu_xxi ))).u()
u_xixi = (poly::sum((lnu_xi_lnu_xi, lnu_xixi))).u()
u_xxx = (poly::sum((poly::times(lnu_x ,lnu_x_lnu_x ), (poly::times(lnu_x , lnu_xx )).C(3), lnu_xxx ))).u()
u_xxxi = (poly::sum((poly::times(lnu_x ,lnu_x_lnu_xi ), (poly::times(lnu_x , lnu_xxi )).C(2),poly::times(lnu_xi, lnu_xx ), lnu_xxxi ))).u()
u_xxixi = (poly::sum((poly::times(lnu_x ,lnu_xi_lnu_xi), (poly::times(lnu_xi, lnu_xxi )).C(2),poly::times(lnu_x, lnu_xixi), lnu_xxixi ))).u()
u_xixixi = (poly::sum((poly::times(lnu_xi,lnu_xi_lnu_xi), (poly::times(lnu_xi, lnu_xixi)).C(3), lnu_xixixi))).u()
u_xxi_xi = u_xxi. dxi()
u_xixi_xi = u_xixi. dxi()
xu_xxi_xi = u_xxi.x().dxi()
if (R>1) {
pointer(class poly scalar) rowvector p; class poly rowvector t; real scalar r
p = &lng_x, &lng_xx, &lng_xxi, &lng_xixi, &lng_xxx, &lng_xxxi, &lng_xxixi, &lng_xixixi
for (i=cols(p); i; i--) {
t = *p[i]
for (r=2; r<=R; r++)
t = t, (t[r-1]).C(1/(r-1)).u() // create sum_r {1/gamma(r)* C * xi^a * sig^b * u^(l+r-1) * (ln u)^d * x^n / y^m}
*p[i] = poly::sum(t)
}
p = &u_x, &u_xx, &u_xxi, &u_xixi, &u_xxx, &u_xxxi, &u_xxixi, &u_xixixi, &u_xxi_xi, &u_xixi_xi, &xu_xxi_xi
for (i=cols(p); i; i--)
*p[i] = (*p[i]).C(1/gamma(R)).u(R-1)
}
}
xlng_xxi = lng_xxi.x()
lng_xxi_xi = lng_xxi.dxi()
lng_xixi_xi = lng_xixi.dxi()
xlng_xxi_xi = xlng_xxi.dxi()
if (GEV) {
lng_x = poly::sum((lng_x ,u_x. C(-1))) // -lny_x + lnu_x - u_x
lng_xx = poly::sum((lng_xx ,u_xx. C(-1))) // -lny_xx + lnu_xx - u_xx
lng_xxi = poly::sum((lng_xxi ,u_xxi. C(-1))) // -lny_xxi + lnu_xxi - u_xxi
xlng_xxi = poly::sum((xlng_xxi,u_xxi.x().C(-1)))
lng_xixi = poly::sum((lng_xixi,u_xixi. C(-1))) // -lny_xixi + lnu_xixi - u_xixi
lng_xxi_xi = poly::sum(( lng_xxi_xi, u_xxi_xi.C(-1)))
xlng_xxi_xi = poly::sum((xlng_xxi_xi,xu_xxi_xi.C(-1)))
lng_xixi_xi = poly::sum((lng_xixi_xi,u_xixi_xi.C(-1)))
lng_xxx = poly::sum((lng_xxx ,u_xxx. C(-1))) // -lny_xxx +lnu_xxx - u_xxx
lng_xxxi = poly::sum((lng_xxxi ,u_xxxi. C(-1))) // -lny_xxxi +lnu_xxxi - u_xxxi
lng_xxixi = poly::sum((lng_xxixi ,u_xxixi. C(-1))) // -lny_xxixi +lnu_xxixi - u_xxixi
lng_xixixi = poly::sum((lng_xixixi,u_xixixi.C(-1))) // -lny_xixixi +lnu_xixixi - u_xixixi
}
lng_xx_xi = lng_xx.dxi()
xlng_x = lng_x.x()
x2lng_xx = lng_xx.x(2)
lng_lnsiglnsig = poly::sum((xlng_x, x2lng_xx))
lng_lnsigxi = xlng_xxi.C(-1)
lng_lnsiglnsig_xi = x2lng_xx.dxi()
lng_lnsigxi_xi = xlng_xxi_xi.C(-1)
lng_lnsiglnsiglnsig = poly::sum((xlng_x.C(-1), x2lng_xx.C(-3), lng_xxx.x(3).C(-1)))
lng_lnsiglnsigxi = poly::sum((xlng_xxi, lng_xxxi.x(2)))
lng_lnsigxixi = lng_xxixi.x().C(-1)
if (GEV) {
xlng_xx = lng_xx.x()
lng_mumu = lng_xx.sig(-2) // lng_xx/sig^2
lng_mulnsig = poly::sum((lng_x, xlng_xx)); lng_mulnsig._sig(-1) // lng_x/sig + xlng_xx/sig
lng_muxi = lng_xxi.sig(-1).C(-1) // -lng_xxi/sig
lng_mumu_xi = lng_xx_xi.sig(-2)
lng_mulnsig_xi = xlng_xx.dxi().sig(-1)
lng_muxi_xi = lng_xxi_xi.sig(-1).C(-1)
lng_mumumu = lng_xxx.sig(-3).C(-1) // -lng_xxx/sig^3
lng_mumulnsig = poly::sum((lng_xx.C(2), lng_xxx.x())); lng_mumulnsig._sig(-2); lng_mumulnsig._C(-1) // -1/sig^2(2lng_xx+xlng_xxx)
lng_mulnsiglnsig = poly::sum((lng_x, xlng_xx.C(3), lng_xxx.x(2))); lng_mulnsiglnsig._C(-1); lng_mulnsiglnsig._sig(-1) // -1/sig*(lng_x+3xlng_xx+xxlng_xxx)
lng_mumuxi = lng_xxxi.sig(-2) // lng_xxxi/sig^2
lng_mulnsigxi = poly::sum((lng_xxi, lng_xxxi.x())); lng_mulnsigxi._sig(-1) // 1/sig*(lng_xxi+xlng_xxxi)
lng_muxixi = lng_xxixi.C(-1); lng_muxixi._sig(-1) // -lng_xxixi/sig
D2lng = lng_mumu , lng_mulnsig , lng_muxi , lng_lnsiglnsig , lng_lnsigxi , lng_xixi
D2lng_xi = lng_mumu_xi , lng_mulnsig_xi , lng_muxi_xi , lng_lnsiglnsig_xi , lng_lnsigxi_xi , lng_xixi_xi
D3lng = lng_mumumu , lng_mumulnsig , lng_mumuxi , lng_mulnsiglnsig , lng_mulnsigxi , lng_muxixi ,
lng_mumulnsig, lng_mulnsiglnsig, lng_mulnsigxi, lng_lnsiglnsiglnsig, lng_lnsiglnsigxi, lng_lnsigxixi,
lng_mumuxi , lng_mulnsigxi , lng_muxixi , lng_lnsiglnsigxi , lng_lnsigxixi , lng_xixixi
} else {
D2lng = lng_lnsiglnsig , lng_lnsigxi , lng_xixi
D2lng_xi = lng_lnsiglnsig_xi , lng_lnsigxi_xi , lng_xixi_xi
D3lng = lng_lnsiglnsiglnsig, lng_lnsiglnsigxi, lng_lnsigxixi,
lng_lnsiglnsigxi , lng_lnsigxixi , lng_xixixi
}
// expand expectations to derivatives w.r.t. coeficients of linear determinants of mu, lnsig, xi
k = X_mu_1+X_lnsig_1+X_xi_1 + cols(X_mu)+cols(X_lnsig)+cols(X_xi)
E_D2lng = J(1,k^2,0); E_D3lng = J(1,k^3, 0)
D = Dmatrix(2+GEV)'
L = Lmatrix(k)'
K = E.E(D2lng)
_E_D2lng = K * D
_E_D3lng = (( (GEV? (J(rows(K),6,0), -2*K[,1],-K[|.,2\.,3|],J(rows(K),3,0)) : J(rows(K),3,0)), E.E(D2lng_xi)) - E.E(D3lng)/2) * I(2+GEV)#D
for (i=rows(_E_D2lng); i; i--) {
Q = blockdiag((X_lnsig[i,],J(1,X_lnsig_1,1)), (X_xi[i,],J(1,X_xi_1,1)))
if (GEV) Q = blockdiag((X_mu[i,],J(1,X_mu_1,1)), Q)
QQ = Q#Q
E_D2lng = E_D2lng + _E_D2lng[i,] * QQ
E_D3lng = E_D3lng + _E_D3lng[i,] * QQ#Q
}
E_D2lng = E_D2lng * L
E_D3lng = E_D3lng * I(k)#L
E_D3lng = colshape(E_D3lng, cols(E_D2lng))'
A = J(k,0,0); for (i=k; i; i--) A = invvech(E_D3lng[,i]), A
return(&E_D2lng, &A)
}
// Cox-Snell bias correction
real rowvector extreme_CS(transmorphic M) {
real scalar i, N, GEV, X_mu_1, X_lnsig_1, X_xi_1, k_mu_covar, k_lnsig_covar, k_xi_covar
real matrix t, A, K, invK, X_mu, X_lnsig, X_xi
pointer(real matrix) rowvector pt
real rowvector b_ML
real colvector subsample, r, ur, lnsig, xi, wt
GEV = moptimize_util_userinfo(M, 1)
t = moptimize_init_userinfo(M, 2); N = rows(t)
r = moptimize_init_userinfo(M, 3)
b_ML = moptimize_result_coefs(M) // ML estimate
X_lnsig_1 = moptimize_init_eq_cons(M, 1+GEV)=="on"
X_xi_1 = moptimize_init_eq_cons(M, 2+GEV)=="on"
X_mu_1 = GEV? moptimize_init_eq_cons(M, 1 )=="on" : 0
t = moptimize_util_eq_indices(M, 1+GEV); k_lnsig_covar = t[2,2]-t[1,2]-X_lnsig_1+1 // number of non-constant covariates
t = moptimize_util_eq_indices(M, 2+GEV); k_xi_covar = t[2,2]-t[1,2]-X_xi_1+1
if (GEV) t = moptimize_util_eq_indices(M, 1 ); k_mu_covar = t[2,2]-t[1,2]-X_mu_1+1
xi = k_xi_covar ? moptimize_util_xb(M, b_ML, 2+GEV) : moptimize_result_eq_coefs(M, 2+GEV)
lnsig = k_lnsig_covar? moptimize_util_xb(M, b_ML, 1+GEV) : moptimize_result_eq_coefs(M, 1+GEV)
X_lnsig = k_lnsig_covar? moptimize_init_eq_indepvars(M, 1+GEV) : J(N,0,0)
X_xi = k_xi_covar? moptimize_init_eq_indepvars(M, 2+GEV) : J(N,0,0)
if (GEV) X_mu = k_mu_covar? moptimize_init_eq_indepvars(M, 1 ) : J(N,0,0)
wt = moptimize_init_weight(M)
if (rows(wt)==1)
wt = cols(X_lnsig) | cols(X_xi)? 1 : N
wt = wt :* (xi:>=-.2)
if (rows(r)==1) {
pt = extreme_CSByR(GEV, r, lnsig, xi, wt, X_mu, X_lnsig, X_xi, X_mu_1, X_lnsig_1, X_xi_1)
K = *pt[1]
A = *pt[2]
} else {
ur = uniqrows(r)
for (i=1; i<=rows(ur); i++) {
subsample = extreme_selectindex(r:==ur[i])
pt = extreme_CSByR(GEV, ur[i], (rows(lnsig)>1? lnsig [subsample,] : lnsig ),
(rows(xi )>1? xi [subsample,] : xi ),
(rows(wt )>1? wt [subsample,] : wt ),
(cols(X_mu )? X_mu [subsample,] : X_mu ),
(cols(X_lnsig)? X_lnsig[subsample,] : X_lnsig),
(cols(X_xi )? X_xi [subsample,] : X_xi ), X_mu_1, X_lnsig_1, X_xi_1)
if (i==1) {
K = *pt[1]
A = *pt[2]
} else {
K = K + *pt[1]
A = A + *pt[2]
}
}
}
invK = invsym(invvech(-K'))
return((invK * A * vec(invK))')
}
// helper function for extreme_BS(): draw and fit one simulated sample
real rowvector extreme_BS_draw(transmorphic scalar M, transmorphic scalar M2, real scalar GEV, real scalar N, real scalar K, real scalar R, real colvector r, real matrix C, real colvector mu, real colvector sig, real colvector xi, real colvector gumbel, real scalar CS) {
real matrix x, x2; real scalar k
do {
x = runiform(N,R)
} while (!all(x)) // generate in range (0,1), not [0,1)
if (GEV) x = R==1? -ln(x) : ln(x)*C
if (rows(xi)==1)
x2 = abs(xi)<1e-10? -ln(x) : (x:^-xi:-1):/xi
else {
x2 = (x:^-xi:-1):/xi
if (rows(gumbel)) x2[gumbel] = -ln(x[gumbel])
}
x2 = x2 * sig; if (GEV) x2 = x2 :+ mu
if (rows(r)>1) // for multiple-order stat models with missing values, copy missingness pattern, moving available order stats to right
for (k=N; k; k--)
if (r[k]1) C = uppertriangle(J(R,R,-1))
if (rows(xi)>1) gumbel = extreme_selectindex(abs(xi):<.001)
bhat = J(reps, K=moptimize_util_eq_indices(M,moptimize_init_eq_n(M))[2,2], .) // width is K, number of params
if (BS2) bhat2 = bhat
for (i=reps; i; i--)
if ( (bhat_BS1 = extreme_BS_draw(M, M2, GEV, N, K, R, r, C, mu, sig, xi, gumbel, CS)) [1] != . )
if (BS2) {
_xi = k_xi_covar ? moptimize_util_xb(M2, bhat_BS1, 2+GEV) : moptimize_result_eq_coefs(M2, 2+GEV)
_lnsig = k_lnsig_covar? moptimize_util_xb(M2, bhat_BS1, 1+GEV) : moptimize_result_eq_coefs(M2, 1+GEV)
if (GEV) _mu = k_mu_covar? moptimize_util_xb(M2, bhat_BS1, 1 ) : moptimize_result_eq_coefs(M2, 1 )
bhat2[i,] = extreme_BS_draw(M2, M2, GEV, N, K, R, r, C, _mu, exp(_lnsig), _xi, J(0,1,0), CS)
if (bhat2[i,1] != .) // only keep 1st bootstrap results if 2nd also attained
bhat[i,] = bhat_BS1
} else
bhat[i,] = bhat_BS1
if (_p != .) {
yp = GEV? -1/ln(1-_p) : 1/_p
xihat = bhat[,2+GEV]
lnzphat = bhat[,1+GEV] :+ ln((yp:^xihat :- 1) :/ xihat)
bhat = bhat, lnzphat
//zphat = exp(meanbhat[1+GEV] + mean(ln(zphat))); if (GEV) zphat = zphat + meanbhat[1]
//meanbhat = meanbhat, zphat
if (BS2) {
xihat2 = bhat2[,2+GEV]
lnzphat2 = bhat2[,1+GEV] :+ ln((yp:^xihat2 :- 1) :/ xihat2)
bhat2 = bhat2, lnzphat2
}
}
meanbhat = mean(bhat); if (BS2) meanbhat2 = mean(bhat2)
V = quadcrossdev(bhat,meanbhat,bhat,meanbhat) / (colnonmissing(bhat)[1] - 1)
return (BS2? 3*meanbhat - 2*b - meanbhat2 : meanbhat - b) // 2nd-order bias estimate is (mean(bhat) - b) - ((mean(bhat2) - mean(bhat)) - (mean(bhat) - b))
}
// Bias correction for EVT models
// reps = 0 for Cox-Snell, otherwise number of replications for parametric bootstrap
// assumes existence of locals bsmall and (if doing BS) Vsmall and b_bs holding names of Stata matrices to take results
void extreme_small(transmorphic M, real scalar reps) {
real rowvector bias, b_bs_mean; real matrix b_bs, V; pragma unset b_bs; pragma unset b_bs_mean
if (reps) {
bias = extreme_BS(M, moptimize_result_coefs(M), reps, 0, 0, V)
st_matrix(st_local("Vsmall"), V)
} else
bias = extreme_CS(M)
st_matrix (st_local("bsmall"), moptimize_result_coefs(M) - bias)
st_matrixcolstripe(st_local("bsmall"), moptimize_result_colstripe(M))
}
/* // helper function for profile plots
// takes moptimize() model, index in full param vector of param to be varied, name of string numlist of evaluation points, offset for displaying numlist values, eq of result to be extracted (eq=0=>log likelihood)
// returns maximized, constrained log likelihoods
real matrix extremeProfile(transmorphic M, real scalar paramindex, string scalar paramname, real colvector numlist, real scalar offset, real scalar ll, real scalar level, real matrix CI) {
real matrix Cns, t; real scalar i, eq
eq=0
(t = J(1, moptimize_util_eq_indices(M, moptimize_init_eq_n(M))[2,2]+1, 0))[paramindex] = 1 // moptimize_util_eq_indices(M, moptimize_init_eq_n(M))[2,2] is just K, number of coefs
Cns = moptimize_init_constraints(M)
Cns = rows(Cns)? Cns \ t : t
moptimize_init_constraints(M, Cns)
moptimize_init_tracelevel(M, "none")
moptimize_init_verbose(M, "off")
moptimize_init_conv_maxiter(M, 100)
t = J(0, 2, 0)
for (i=1; i<=rows(numlist); i++) {
Cns[rows(Cns),cols(Cns)] = numlist[i]
moptimize_init_constraints(M, Cns)
if (_moptimize(M))
printf("Unable to maximize profile likelihood for %s = %f.\n", paramname, numlist[i]+offset)
else
t = t \ (numlist[i], (eq? moptimize_result_eq_coefs(M,eq) : moptimize_result_value(M)))
}
return(t)
}
*/
real scalar extreme_r0_to_ll(transmorphic M, real matrix Cns, real scalar r0) {
Cns[rows(Cns),cols(Cns)] = r0
moptimize_init_constraints(M, Cns)
if (_moptimize(M))
return (.)
return (moptimize_result_value(M))
}
real scalar extremeSearch(transmorphic M, real matrix Cns, real scalar alpha, real scalar p_lo, real scalar lo, real scalar p_hi, real scalar hi) {
real scalar mid, _p
mid = (alpha-p_lo)/(p_hi-p_lo)*(hi-lo) + lo
if (mreldif(lo,mid)<1e-6 | mreldif(hi,mid)<1e-6)
return (mid)
if ( ((_p = extreme_r0_to_ll(M, Cns, mid)) < alpha) == (p_lo < alpha) )
return(extremeSearch(M, Cns, alpha, _p, mid, p_hi, hi))
return(extremeSearch(M, Cns, alpha, p_lo, lo, _p, mid))
}
// derive profile-based CI. Temporarily modifies passed moptimize() estimation object by adding a constraint.
// returns Y coords of profile plot. Return argument CI holds bounds of possibly disjointed confidence set, one range per row.
real matrix extremeProfile(transmorphic M, real scalar paramindex, string scalar paramname, real colvector plotX, real scalar offset, real scalar ll, real scalar level, real matrix CI) {
real scalar alpha, i, j, lo, hi, se, z_alpha, b, missing; real colvector plotY; real matrix _Cns, Cns, t
moptimize_init_tracelevel(M, "none")
moptimize_init_verbose(M, "off")
moptimize_init_conv_warning(M, "off")
moptimize_init_search(M, "off")
moptimize_init_conv_maxiter(M, 100)
(t = J(1, cols(moptimize_init_coefs(M))+1, 0))[paramindex] = 1 // make new constraint with coef of 1 on constrained parameter; moptimize_util_eq_indices(M, moptimize_init_eq_n(M))[2,2] is just K, number of coefs
Cns = _Cns = moptimize_init_constraints(M)
Cns = rows(Cns)? Cns \ t : t
alpha = ll - invchi2(1, level*.01) * .5
if (plotX == .) {
if (_moptimize(M)) { // fit unconstrained model once to perform Wald test for starting points
printf("Unable to maximize likelihood for %s.\n", strlower(paramname))
return (.)
}
b = moptimize_result_coefs(M)[paramindex]
se = sqrt(moptimize_result_V(M)[paramindex,paramindex])
z_alpha = abs(invnormal((1-level*.01)*.5))
lo = b - z_alpha * se
hi = b + z_alpha * se
for (i=10; i & -extreme_r0_to_ll(M, Cns, lo)<-alpha; i--) // try at most 10 times to bracket confidence set by doubling on high and low sides in turn
lo = 2 * lo - hi
for (i=10; i & -extreme_r0_to_ll(M, Cns, hi) < -alpha; i--)
hi = 2 * hi - lo
plotX = rangen(lo, hi, 25)
}
plotY = J(rows(plotX), 1, .)
for (i=rows(plotX); i; i--)
if (missing = (. == (plotY[i] = extreme_r0_to_ll(M, Cns, plotX[i]))))
printf("Unable to maximize profile likelihood for %s = %f.\n", strlower(paramname), plotX[i]+offset)
if (level<100) {
if (missing) {
t = plotY:<.
plotX = select(plotX, t)
plotY = select(plotY, t)
}
if (rows(plotX)) {
CI = plotY :> alpha; CI = CI[|2\.|] - CI[|.\rows(plotX)-1|]
lo = extreme_selectindex(CI:== 1)
hi = extreme_selectindex(CI:==-1)
if (rows(lo)==0 & rows(hi)==0)
CI = . , .
else {
if (rows(lo)==0) lo = .
else if (rows(hi)==0) hi = .
else {
if ( lo[1 ] > hi[1 ]) lo = . \ lo // non-rejection ranges that are not bounded within grid
if (-lo[rows(lo)] < -hi[rows(hi)]) hi = hi \ .
}
CI = lo, hi
for (i=rows(lo); i; i--)
for (j=2; j; j--)
if (CI[i,j]<.)
CI[i,j] = extremeSearch(M, Cns, alpha, plotY[CI[i,j]], plotX[CI[i,j]], plotY[CI[i,j]+1], plotX[CI[i,j]+1])
}
}
}
moptimize_init_constraints(M, _Cns)
return ((plotX,plotY))
}
// helper function for MRL plot--fast computation of mean exceedence at each data point, and standard error thereof
// assumes X is sorted ascending
/*real matrix extreme_mrl(string scalar Xname, string scalar wtname) {
real colvector X, wt, mu, se; real scalar N, i, t, dmu
st_view(X, ., Xname)
st_view(wt, ., wtname)
wt = wt / sum(wt)
mu = se = J(N=rows(X), 1, 0)
mu[N] = X[N]
for (i=N-1; i; i--) {
t = X[i] - mu[i+1]
mu[i] = mu[i+1] + (dmu = (X[i]-mu[i+1])/(N-i+1))
se[i] = se[i+1] + dmu*dmu*(N-i) + t*t
}
return (mu-X, sqrt(se):/(N-1::0))
}*/
// helper function for MRL plot--fast computation of mean exceedence at given data point, and CI thereof
// assumes X is sorted ascending
real rowvector extreme_mrl(string scalar Xname, string scalar wtname, string scalar wttype, real scalar min, real scalar level) {
real colvector _X, X, wt, t, select; real scalar N, mu, se, sumwt, halfci
st_view(_X, ., Xname)
select = _X:>=min
st_view(X , extreme_selectindex(select), Xname)
st_view(wt, extreme_selectindex(select), wtname)
if (cols(wt)) {
N = wttype=="fweight"? colsum(wt) : rows(X)
sumwt = colsum(wt)
mu = (X ' wt) / sumwt
t = X :- mu
se = sqrt((t:*t) ' wt / sumwt / (N - 1))
} else {
N = rows(X)
mu = colsum(X) / N
t = X :- mu
se = sqrt(colsum(t:*t) / N / (N-1))
}
mu = mu - min
halfci = se * invttail(N-1, .5 -level/200)
return (mu, mu-halfci, mu+halfci)
}
mata mlib create lextreme, dir("`c(sysdir_plus)'l") replace
mata mlib add lextreme *(), dir("`c(sysdir_plus)'l")
mata mlib index
end