*! version 1.2.0 17 May 2017
help extreme

Title

Extreme value theory models: fitting and graphical diagnostics

Syntax

extreme model varlist [if] [in] [weight] [, estimation_options]

- or, usually after the above -

extreme plot [varname] [if] [in] [weight] [, plot_options]

The varlist in an extreme model command line should have one entry unless fitting a model with multiple order statistics. The varname in an extreme plot command
is only used for the mrl plot type. All varnames and varlists may include factor variable and time-series operators; see fvvarlist.

model

gpd Generalized Pareto distribution (GPD)

gev Generalized Extreme Value distribution (GEV) with one or more order statistics

estimation_options

Model
threshold(varname | numlist) Set minimum value(s) for estimation sample, relative to which excedences are computed; GPD only

muvars([varlist] [, noconstant]) Specify linear determinants (covariates) for mu (location) parameter; GEV only

sigvars([varlist] [, noconstant]) Specify linear determinants for log sigma (scale) parameter

xivars([varlist] [, noconstant]) Specify linear determinants for xi (shape) parameter

gumbel Impose xi=0: exponential distribution for GPD, Gumbel distribution for GEV

constraints(constraints) Impose other linear constraints on parameters

Estimation
init(rowvector) Override default starting point for maximization; see help ml.

ml_opts Standard ml maximize options to control Maximum Likelihood search; seldom needed

Corrections to classical estimates
vce(vcetype) vcetype may be robust; cluster clustvar; opg; oim; bootstrap [, bootstrap_options]; or (equivalently) bs [,
bootstrap_options]; see vce_option

small({cs | bs [, reps(#)]}) cs: Analytical (Cox-Snell) bias correction to coefficient estimates

bs: parametric-bootstrap-based bias correction and computation of VCE

Display
level(#) Set confidence level; default is 95

quietly Suppress most output other than final results

display_opts Standard ml display options to control results display; seldom needed

Graphing
graph_opts Most graph twoway options for formatting plots; relevant for GPD estimation with multiple thresholds specified, which triggers graphing

plot_options

mrl[(numlist)] Mean residual life diagnostic plot for varname; may be requested before estimation

pp Probability-probability diagnostic plot

qq Quantile-quantile diagnostic plot

return Return level-return period diagnostic plot

density Kernel density diagnostic plot

xiprofile([numlist]) Profile likelihood plot for xi at specified values

retprofile([numlist], period(#)) Profile likelihood plot for return level, for specified return period

rateprofile([numlist], retlevel(#)) Profile likelihood plot for return rate, for specified return level

level(#) Set confidence level for plotted confidence intervals; default is 95

name(stub) Specify prefix for storage names of plots produced; default is extreme

graph_opts Most graph twoway options for formatting plots

kdens_opts When requesting the density plot, standard twoway kdensity options may be included

extreme shares features of all estimation commands; see help estcom

The syntax of predict following extreme is

predict [type] newvarname [if] [in] [, statistic equation(#eqno|eqname)]

where statistic is xb, pr, cdf, ccdf, or invccdf(# | varname | exp)

Donate?

Has extreme improved your career or your marriage?
Consider giving back through a donation to support the work of its self-employed author, David Roodman.

Description

extreme estimates the most-used models in univariate extreme value theory, using Maximum Likelihood. It is strongly influenced by the textbook of Coles (2001)
and the associated ismev package for R, maintained by Eric Gilleland. In addition to fitting models,
extreme offers several diagnostic
and profile plots. As a result, and as a measure of its scope, extreme can run most of the numerical and graphical examples in Coles (2001). (See examples below.)

Specifically, extreme fits:

* The generalized Pareto distribution (GPD), which is appropriate for modeling exceedances of a threshold, i.e., right tails of distributions. The GPD cumulative distribution function (cdf)
takes the form G(x)=1-(1+xi*(x-mu)/sigma)^(-1/xi), where mu is a pre-specified threshold.

* The generalized extreme value distribution (GEV), which is appropriate for modeling block maxima. Its cdf is G(x)=exp(-(1+xi*(x-mu)/sigma)^(-1/xi)).

* The extension of the GEV for multiple order statistics for blocks, such as the highest, second-highest, and third-highest daily rainfall totals in a
year. Its cdf does not have a closed form, but its probability density function (pdf) does (Coles 2001, eq 3.15). The number of order statistics available may vary by observation.

The GPD model has two parameters to be estimated: a scale parameter, sigma, and a shape parameter, xi. extreme enters sigma in logarthmic form, so that the
core parameter displayed in estimation results is lnsigma. (Normally the implied value for sigma is also displayed, at the bottom.) Fitting
the GPD also requires specifying a threshold, relative to which exceedances are computed. Observations at or below the threshold are excluded.

The GEV model too contains scale and shape parameters. In addition, in place of the pre-specified threshold, there is a third parameter for location, mu.

Multiple order statistics are specified by including multiple variables in the varlist after extreme gev, in descending order. For each observation,
extreme will model as many order statistics as available data permit. E.g., if you list 10 variables and just variable 5 is missing for some observation,
then extreme will model the top 4 statistics for that observation.

extreme allows for non-stationary models: all core parameters may depend linearly on variables or factor variables specified with the muvars(), sigvars(), xivars()
options.

extreme's major novelty is the ability to compute the Cox-Snell small-sample bias correction for all models fit, via the small() option. The correction for the GPD is
derived in Giles, Feng, and Godwin (2015). The correction for the GEV, including for multiple order statistics, is new. Its efficacy has been demonstrated
in simulations, and a write-up is pending. The correction has been extended to the non-stationary case too. Maximum Likelihood is always biased in finite samples, and the bias can be
significant in the small samples often used in extreme value analysis. (On analytical bias correction generally, see Cordeiro and Cribaro-Neto (2014, ch. 4).)

By default, the variance-covaraince matrix
from the classical Maximum Likelihood fit and, thus, the standard errors are retained after this correction to the parameter estimates. Simulations suggests that the correction hardly affects
the standard errors. But you can at once bypass and check this assertion by bootstrapping the standard errors of the corrected estimator, e.g., by adding a vce(bs) option.

However, extreme also offers bias-correction via a parametric bootstrap. This has the advantage of working for values of xi for which ML is consistent, but for which the
Cox-Snell correction is not defined: formally between -1 and -1/3, in practice between -1 and about -0.2. Also, it works for models with constraints, including the gumbel
constraint. When bootstrapping
the bias correction, extreme also automatically bootstraps the variance-covariance matrix of the estimate.

When bootstrapping the VCE or bias correction, the precise results depend on the initial state of Stata's pseudorandom number generator. See set seed. extreme displays this state
and stores it in the return value e(seed).

As an estimation command built on Stata's ml, extreme accepts observation weights, linear constraints on parameters, vce() options, and options
for controlling optimization, all according to standard ml syntax.

When fitting the GPD, specifying multiple thresholds in the threshold() option causes extreme to fit the model for each threshold, adjusting the sample accordingly, and then
graph the estimates of lnsig and xi against the
thresholds. In addition, all of the extreme plot subcommands of course generate plots. In all cases, extreme makes minimal effort to format and label
graphs. This reduces interference with any graph twoway options that you include in the command line, to change the scheme, titles, etc.

You can request multiple plots in one extreme plot command line, as in extreme plot, qq pp density xiprofile. Most plot types do not affect the stored numerical estimation
results. However, requesting xiprofile(), retprofile(), or rateprofile() will cause extreme to save matrices named e(xiprofileCI) and e(retprofileCI) holding the confidence
sets derived through graphing the profile likelihood.

Anyone working with these models is likely to be familiar with Coles (2001). The easiest way to learn extreme is probably to study the replications below
of many examples in that book.

Estmation options

threshold(varname | numlist) is for the GPD model only. It sets the minimum value for inclusion in the sample. For fitting,
extreme subtracts the threshold from all observations above it. Observations less than or equal to it are excluded. The threshold defaults to the minimum value in the sample.
If a varname is specified, then the threshold may vary by observation.
If multiple thresholds are specified, extreme graphs the estimated lnsig and xi parameters against the chosen thresholds. If exactly two thresholds are specified,
extreme interprets them as minimum and maximum values with 100 steps in between.

muvars([varlist] [, noconstant]), sigvars([varlist] [, noconstant]),
xivars([varlist] [, noconstant]) specify covariates for non-stationary models. They allow mu, log sigma (not sigma), and xi
to depend linearly on covariates. Factor variables are allowed.

gumbel imposes the constraint xi=0. In the context of the GPD model, this specifies the exponential model. In the context of the GEV, this specifies
the Gumbel model.

constraints(constraints) imposes constraints according to standard Stata syntax.

level(#) specifies the confidence level, in percent, for confidence intervals; see help level. The default is 95. It affects display of regression results
and plots of confidence intervals.

quietly suppresses most output, but not final results.

init(vector) passes a row vector of user-chosen starting values for the model fit, in the manner of the ml init, copy
command. The vector must contain exactly one element for each parameter extreme will estimate, and in the same order as in the output (excluding the derived sigma results at bottom). The default initial coefficient for any non-constant
covariate is 0. As for constants--including all parameters in stationary models--the default initial value for the xi equation is also 0. The defaults for lnsigma and, in the GEV model, mu are what the true values would be if xi=0 and
the observed first and second moments of the modeled variable exactly equal the true ones. That is, fixing xi=0, if x ~ GPD, then Var[x]=sigma^2; and if x ~ GEV,
E[x] = mu + sigma * gamma, Var[x] = sigma^2*pi^2/6, where gamma is the Euler-Mascheroni constant. Solving these equations for mu and sigma and substituting empirical values
for E[x] and Var[x] yields the default initial values.

ml_opts: extreme accepts the following standard ml options to control the search for the maximum: trace
nolog
gradient
hessian
showstep
technique(algorithm_specs)
iterate(#)
tolerance(#)
ltolerance(#)
gtolerance(#)
nrtolerance(#)
nonrtolerance
shownrtolerance
difficult

vce(oim | opg | robust | cluster varname | {bootstrap | bs} [, bootstrap_options]) governs the estimation of the
variance-covariance matrix of the parameter estimates (VCE). See vce_option. If a bootstrapped VCE is chosen, the bootstrap is non-parametric, as is standard in Stata, meaning that
the data are randomly resampled for each replication. This contrasts with the parametric bootstrap deployed by small(bs), discussed next.

small({cs | bs [, reps(#)]}) requests bias correction. ML has o(1/N) bias, which can be significant for small samples. small(cs) requests the
analytical Cox-Snell (1968) bias correction; this is available for xi>-0.2, so if xi<=-.2 no correction is made. The CS correction for
the GPD is derived in Giles, Feng, and Godwin (2015). The correction for the GEV, including for multiple order statistics, is from Roodman (2017).

Alternatively, small(bs) requests a simulation-based bias correction. A parametric bootstrap
is performed--"parametric" meaning that instead of resampling the data as with vce(bs), simulations of the modeled variable are repeatedly drawn from the
ML-estimated distribution (Horwitz 2001, section 3.1). small(bs) also automatically bootstraps the VCE during this process. The optional reps(#) suboption controls the number of replications; the default is
50. As an example, suppose ML estimates xi=0.27 on a
small sample. It is
then found that in simulated data sets of the same size with xi=0.27, the ML estimate averages 0.29, for an empirical bias of +0.02. That suggests that the bias embedded in our
original ML estimate is also about 0.02, so we correct it to 0.25.

Combining small(bs) with vce(bs) tells extreme to nonparametrically bootstrap the VCE of the parametric-bootstrap-corrected coefficients (and thus to set aside
the parametrically bootsrapped VCE computed anyway). This is feasible but slow.

display_opts: extreme accepts the following standard ml display options: noheader
nofootnote
first
neq(#)
showeqns
plus
nocnsreport
noomitted
vsquish
noemptycells
baselevels
allbaselevels
cformat(%fmt)
pformat(%fmt)
sformat(%fmt)
nolstretch
coeflegend

graph_opts: Most graph twoway options for formatting plots can be added to extreme gpd command lines
when multiple thresholds are provided, in order to override formatting defaults.

Plot options

An extreme plot command line may include multiple plot types.

mrl[(numlist)] requests a mean residual life diagnostic plot, which
graphs the average amount by which observations above a threshold
exceed it, as a function of the threshold. If the tail of a distribution is well-approximated by the GPD
above some threshold, then the MRL plot should be linear above that point
(Coles 2001, section 4.3.1). Among extreme plot types, mrl is anomolous in being available
before estimation and, by the same token, requiring a
varname before the command line's main comma. The numlist indicates at which
thresholds the mean residual should be computed. If the numlist
is omitted, all observed values of the variable
are used, which can be time-consuming in large data sets. If exactly two values are provided,
these are taken as the minimum and maximum with 100
steps in between; missing values for either bound are interpreted as the variable's minimum or
maximum. Otherwise, the mean residual is computed only at the listed values.

pp plots modeled versus empirical cumulative distribution functions. If the model is
a good fit, the plotted points should track closely to a 45-degree
line. See [R] diagnostic plots and Coles (2001, section 3.3.5). The "empirical cumulative
distribution function" is essentially the (weighted) share of the observations smaller or equal to
a given observation. E.g, if there are 99 observations, the values of the empirical cdf are 0.01,
0.02...0.99. (This formulation avoids working with the empirical cdf to 0 or 1, which could
correspond to infinite values of the modeled variable.) The modeled cdf for each observation is
computed using the estimated parameters.

qq plots modeled versus empirical cumulative quantiles of the distribution of the modeled
variable. Again, if the model fits well, the plotted points should follow a 45-degree
line. See [R] diagnostic plots and Coles (2001, section 3.3.5). The empirical quantiles are
just the values of the modeled variable, in order. The modeled quantiles result from
applying the inverse cumulative distribution function to "empirical cdf" values listed above, using
the estimated parameters. For non-stationary models, the data are
pre-transformed to the Gumbel distribution, as described in Coles (2001, section 6.2.3).

return requests a graph of the return level as a function of the return period, as
predicted by the fitted model, and compares that to the pattern of extremes
in the data. See Coles (2001, section 3.3.5). E.g., if the model fit implies that a category 4 or
higher hurricane occurs every 20 years, then one should appear in the data about
that often. Available for stationary models only.

density produces a plot comparing the modeled distribution to the actual distribution, the latter respresented by
superimposed dot and kernel density plots. Available for stationary models only.

xiprofile([numlist]) requests a likelihood profile plot for the shape parameter xi. The model is
repeatedly estimated with xi each time constrained
to a different value. The
maximized (profile) log likelihood is plotted against xi. This can provide a more reliable
basis for constructing confidence intervals. The coverage of the confidence interval
is controlled by the level() option. The derived confidence interval is stored in
the matrix e(xiprofileCI). The optional numlist is parsed as with the
mrl option, except that missing values are not allowed. If the numlist is omitted,
extreme will automatically choose endpoints for the search. Available only for models
that are stationary in xi.

retprofile([numlist], period(#)), analogous to
the previous option, requests a likelihood profile plot for the return level for a given
return period. The model is
repeatedly estimated with the return level each time constrained to a different value. The
maximized (profile) log likelihood is plotted against the return level. The derived confidence
interval is stored in the matrix e(retprofileCI). The required
suboption period(#) specifies the return period. Available for stationary models only.

rateprofile([numlist], retlevel(#)) requests a
likelihood profile plot for the return rate (reciprocal of return period), for a given return level. The model is
repeatedly estimated with the return rate each time constrained to a different value. The
maximized (profile) log likelihood is plotted against the return level. The derived confidence
interval is stored in the matrix e(rateprofileCI). The required
suboption retlevel(#) specifies the return magnitude. Available for stationary models only.

graph_opts: Most graph twoway options for formatting plots can be added to extreme plot command lines
to override plot formatting defaults, from schemes, to fonts, to scaling, to labelling.

kdens_opts: When running the density plot, standard twoway kdensity options may be included to tune the kernel regressions.

predict syntax

Options for predict after extreme are:

xb The default: linear prediction of mu, lnsig, or xi

equation(#eqno|eqname) Equation for linear prediction; irrelavant for other statistics

pr Probability density for each observation

cdf Cumulative probability density

ccdf Complementary cumulative probability density (1-cdf, return rate)

invccdf(#|varname|exp) Inverse complementary cumulative probability density (return level) for return rate specified by a constant, variable, or expression

Stored results

extreme stores the following in e():

Scalars

e(rank) rank of e(V)

e(ic) number of iterations

e(k) number of parameters

e(k_eq) number of equations (2 for GPD, 3 for GEV)

e(k_dv) number of dependent variables (same)

e(converged) 1 if converged, 0 otherwise

e(rc) return code

e(k_eq_model) number of equations in overall model test

e(df_m) model degrees of freedom

e(chi2) statistic from overall model fit test

e(p) associated p value

e(ll) log likelihood

e(Nthresh) number of thresholds in thresh() option (for GPD)

e(Ndepvar) number of dependent variables (>1 for GEV with multiple order statistics)

e(zeta) fraction of potential sample exceeding threshold (for GPD)

e(stationary) 0 if model has covariates (in muvars(), sigvars(), xivars() options), 1 otherwise

Macros

e(cmd) "extreme"

e(cmdline) command as typed

e(predict) program used to implement predict

e(muvars) content of muvars() option, if any

e(sigvars) content of sigvars() option, if any

e(xivars) content of xivars() option, if any

e(diparmopts) options passed to ml display when displaying results

e(gumbel) "gumbel" if that option used

e(model) model type (cmd:gpd} or gev

e(depvar) modeled variable

e(seed) initial pseudorandom number generator state, if bootstrap estimation

e(chi2type) "Wald"

e(opt) "moptimize"

e(vce) vce() type

e(title) title for results

e(user) program to compute likelihood

e(ml_method) ml optimization method (lf2)

e(technique) search technique

e(which) max

e(properties) b V

Matrices

e(b) parameter vector

e(V) variance-covariance matrix of the parameter estimates

e(ilog) iteration log

e(gradient) gradient vector

e(xiprofileCI) profile-based confidence interval for xi (after extreme plot, xiprofile())

e(xiprofileplot) matrix containing coordinates of graphed profile for xi

e(retprofileCI) profile-based confidence interval for return level (after extreme plot, retprofile())

e(retprofileplot) matrix containing coordinates of graphed profile for return level

e(rateprofileCI) profile-based confidence interval for return rate (after extreme plot, rateprofile())

e(rateprofileplot) matrix containing coordinates of graphed profile for return rate

Function

e(sample) marks estimation sample

Examples

All examples are clickable.

* Examples from Coles (2001), chapter 3

. use http://fmwww.bc.edu/repec/bocode/p/portpirie

. extreme gev SeaLevel

. extreme plot, xiprofile(-.3 .3) retprofile(4.15 4.6, period(10)) // Figures 3.2, 3.3

. extreme plot, retprofile(4.4 6, period(100)) pp qq return density // Figures 3.4, 3.5

. est store GEV // save results for this unrestricted model

. extreme gev SeaLevel, gumbel // impose Gumbel constraint (xi = 0)

. lrtest GEV // likelihood ratio test of Gumbel vs. unrestricted model: accepts restriction

. extreme plot, pp qq return density // Figure 3.6

. use http://fmwww.bc.edu/repec/bocode/g/glass

. gen negBreakingStrength = -BreakingStrength

. extreme gev negBreakingStrength

. extreme plot, pp qq return density // Figure 3.7

. use http://fmwww.bc.edu/repec/bocode/v/venice

. extreme gev SeaLevel1 // Table 3.1

. extreme gev SeaLevel1-SeaLevel5 // Table 3.1

. extreme plot, pp qq return density // Figure 3.9

. extreme gev SeaLevel1-SeaLevel10 // Table 3.1: use 10 order stats, except only 6 availble for 1935

. forvalues r=2/10 { (Click here to run loop)

. extreme gev SeaLevel1-SeaLevel`r'

. extreme plot, return name(v`r') nodraw title(r=`r')

. }

. graph combine v2Return v3Return v4Return v5Return v6Return v7Return v8Return v9Return v10Return // Figure 3.8

* Examples from Coles (2001), chapter 4

. use http://fmwww.bc.edu/repec/bocode/r/rain

. extreme plot Rainfall, mrl(. .) // Figure 4.1

. extreme gpd Rainfall, thresh(0(2.5)50) // Figure 4.2

. extreme gpd Rainfall, thresh(30)

. extreme plot, pp qq return density // Figure 4.3

. nlcom Return100: `e(threshold)' + exp([lnsig]_cons)/[xi]_cons * ((365*100 * `e(zeta)')^[xi]_cons - 1) // conf. interval for 100-year return level

. extreme plot, xiprofile(-.1 .6) retprofile(79 220, period(36500)) // Figures 4.4, 4.5

. use http://fmwww.bc.edu/repec/bocode/d/dowjones

. gen change = 100*ln(Index/Index[_n-1])

. extreme plot change, mrl(-7 5) // Figure 4.6

. extreme gpd change, thresh(2)

. extreme plot, pp qq return density // Figure 4.7

* Examples from Coles (2001), chapter 6

. use http://fmwww.bc.edu/repec/bocode/p/portpirie

. extreme gev SeaLevel

. est store Stationary

. extreme gev SeaLevel, muvar(Year) // add covariate Year to location equation

. lrtest Stationary // LR test comparing it to stationary model

. use http://fmwww.bc.edu/repec/bocode/f/fremantle

. replace Year = Year - 1896 // code 1897 as 1

. extreme gev SeaLevel // stationary model

. est store Stationary

. extreme gev SeaLevel, muvar(Year) // add covariate Year to location equation

. predict muhat

. scatter SeaLevel Year || line muhat Year // Figure 6.1

. lrtest Stationary // LR test comparing it to stationary model

. est store muYear

. extreme plot, pp qq // Figure 6.2

. extreme gev SeaLevel, muvar(SOI) // add covariate Southern Oscillation Index to location equation

. extreme gev SeaLevel, muvar(Year SOI) // add Year back too

. lrtest muYear // LR test comparing this to Year-only model

. use http://fmwww.bc.edu/repec/bocode/v/venice

. extreme gev SeaLevel*, muvar(Year)

. predict muhat

. scatter SeaLevel* Year, pstyle(p1...) msize(vsmall...) || line muhat Year, legend(off) // Figure 6.5

. use http://fmwww.bc.edu/repec/bocode/r/rain

. gen time = _n

. extreme gpd Rainfall, thresh(30) sigvar(time) // model log sigma = beta_0 + beta_1*time

. extreme plot, pp qq // Figure 6.7

. use http://fmwww.bc.edu/repec/bocode/w/wooster

. gen negTemperature = -Temperature

. gen byte Season = mod(floor((Day+30)/(365.25/4)),4)+1 // seasons definition?

. mat threshmat = -10, -25, -50, -30 // Coles seasonal thresholds

. gen threshvar = threshmat[1,Season]

. bysort Season: extreme gpd negTemperature, thresh(threshvar) // separate model for each season

. extreme gpd negTemperature, thresh(threshvar) sigvar(ibn.Season, nocons) // combined model with single shape parameter

* Other examples

. use http://fmwww.bc.edu/repec/bocode/p/portpirie

. extreme gev SeaLevel

. extreme plot, xiprofile(-.3 .3) retprofile(4.15 4.6, period(10)) rateprofile(.05 .3, retlevel(4.2))

. mat list e(xiprofileCI) // programmatically extract profile-based confidence interval for xi

. mat list e(retprofileCI) // programmatically extract profile-based confidence interval for return level

. mat list e(rateprofileCI) // programmatically extract profile-based confidence interval for return rate

. extreme plot, qq scheme(s1rcolor) aspect(1) title(Quantile-quantile plot) xtitle(Modeled sea level quantile) // Q-Q plot with altered appearance

. extreme gev SeaLevel, muvar(Year) small(cs) // Cox-Snell bias correction for non-stationary GEV model

. extreme gev SeaLevel, muvar(Year) small(cs) vce(bs, reps(100)) // Same, with bootstrapped standard errors

. extreme gev SeaLevel, muvar(Year) small(bs, reps(200)) // Parametric-bootstrap-based bias correction and standard errors

. predict mypr, pr // Predict probability densities

. predict myretlevel, invccdf(1/100) // Predict 100-year SeaLevel value, as modeled function of Year

References

Coles, S. 2001. An introduction to statistical modeling of extreme values. Springer.

Cordeiro, G., and F. Cribaro-Neto. 2014. Analytical and bootstrap bias correctlions. In Introduction to Bartlett Correction and Bias Reduction. Springer.

Cox, D.R, and E.J. Snell. 1968. A general definition of residuals. Journal of the Royal Statistical Society. Series B.

Giles, D.A., H. Feng, and R.T. Godwin. 2015. Bias-corrected maximum likelihood estimation of the parameters of the generalized Pareto
distribution. Communications in Statistics - Theory and Methods.

Horwitz, J.L. 2001. The bootstrap. Handbook of Econometrics. Elsevier Science.

Roodman, D. 2017. Bias and size corrections in extreme value modeling. Communications in Statistics - Theory and Methods. DOI: 10.1080/03610926.2017.1353630.

Author

David Roodman

david@davidroodman.com