Stata package help file
Note
This page simply reproduces, for convenience, the output of typing help edm into the Stata console.
The command edm implements a series of tools that can be used for empirical dynamic modeling in Stata. The core algorithm is written in C++ (with a Mata backup) to achieve reasonable execution speed. The command keyword is edm, and should be immediately followed by a subcommand such as explore or xmap. A dataset must be declared as time-series or panel data by the tsset or xtset command prior to using the edm command, and time-series operators including l., f., d., and s. can be used (the last for seasonal differencing).
Syntax
The explore subcommand follows the syntax below and supports one variable for exploration using simplex projection or S-mapping.
edm explore variable [if exp], [e(numlist ascending >=2)]
[tau(integer 1)] [theta(numlist ascending)] [k(integer 0)] [ALGorithm(string)] [REPlicate(integer 0)]
[seed(integer 0)] [full] [RANDomize] [PREDICTionsave(name)] [COPREDICTionsave(name)] [copredictvar(string)]
[CROSSfold(integer 0)] [CI(integer 0)] [EXTRAembed(string)] [ALLOWMISSing] [MISSINGdistance(real 0)]
[dt] [reldt] [DTWeight(real 0)] [DTSave(name)] [DETails] [reportrawe] [strict] [Predictionhorizon(string)]
[dot(integer 1)] [mata] [gpu] [nthreads(integer 0)] [savemanifold(name)] [idw(real 0)]
[lowmemory]
The second subcommand xmap performs convergent cross-mapping (CCM). The subcommand follows the syntax below and requires two variables to follow immediately after xmap. It shares many of the same options with the explore subcommand although there are some differences given the different purpose of the analysis.
edm xmap variables [if exp], [e(integer 2)] [tau(integer 1)] [theta(real 1)]
[Library(numlist)] [RANDomize] [k(integer 0)] [ALGorithm(string)] [REPlicate(integer 0)] [strict]
[DIrection(string)] [seed(integer 0)] [PREDICTionsave(name)] [COPREDICTionsave(name)] [copredictvar(string)]
[CI(integer 0)] [EXTRAembed(string)] [ALLOWMISSing] [MISSINGdistance(real 0)] [dt] [reldt]
[DTWeight(real 0)] [DTSave(name)] [oneway] [DETails] [SAVEsmap(string)] [Predictionhorizon(string)]
[dot(integer 1)] [mata] [gpu] [nthreads(integer 0)] [savemanifold(name)] [idw(real 0)]
[lowmemory]
The third subcommand update updates the plugin to its latest version
The fourth subcommand version displays the current version number
Options
Options for explore and xmap subcommands
e(numlist ascending): This option specifies the number of dimensions E used for the main variable in the manifold reconstruction. If a list of numbers is provided, the command will compute results for all numbers specified. The xmap subcommand only supports a single integer as the option whereas the explore subcommand supports the option as a numlist. The default value for E is 2, but in theory E can range from 2 to almost half of the total sample size. The actual E used in the estimation may be different if additional variables are incorporated. A error message is provided if the specified value is out of range. Missing data will limit the maximum E under the default deletion method.
tau(integer): The tau (or τ) option allows researchers to specify the ‘time delay’, which essentially sorts the data by the multiple τ. This is done by specifying lagged embeddings that take the form: t,t-1τ,…,t-(E-1)τ, where the default is tau(1) (i.e., typical lags). However, if tau(2) is set then every-other t is used to reconstruct the attractor and make predictions—this does not halve the observed sample size because both odd and even t would be used to construct the set of embedding vectors for analysis. This option is helpful when data are oversampled (i.e., spaced too closely in time) and therefore very little new information about a dynamic system is added at each occasion. However, the tau() setting is also useful if different dynamics occur at different times scales, and can be chosen to reflect a researcher’s theory-driven interest in a specific time-scale (e.g., daily instead of hourly). Researchers can evaluate whether τ>1 is required by checking for large autocorrelations in the observed data (e.g., using Stata’s corrgram function). Of course, such a linear measure of association may not work well in nonlinear systems and thus researchers can also check performance by examining ρ and MAE at different values of τ.
theta(numlist ascending): Theta (or θ) is the distance weighting parameter for the
local neighbours in the manifold. It is used to detect the nonlinearity of the system in the explore
subcommand for S-mapping. Of course, as noted above, for simplex projection and CCM a weight of
theta(1) is applied to neighbours based on their distance, which is reflected in the fact that the
default value of θ is 1. However, this can be altered even for simplex projection or CCM (two cases
that we do not cover here). Particularly, values for S-mapping to test for improved predictions as
they become more local may include the following command: theta(0 .00001 .0001 .001 .005 .01 .05 .1 .5 1 1.5 2 3 4 6 8 10).
k(integer): This option specifies the number of neighbours used for prediction. When
set to 1, only the nearest neighbour is used, but as k increases the next-closest nearest neighbours
are included for making predictions. In the case that k is set 0, the number of neighbours used is
calculated automatically (typically as k = E + 1 to form a simplex around a target), which is the
default value. When k < 0 (e.g., k(-1)), all possible points in the prediction set are used (i.e.,
all points in the library are used to reconstruct the manifold and predict target vectors). This
latter setting is useful and typically recommended for S-mapping because it allows all points in the
library to be used for predictions with the weightings in theta. However, with large datasets this
may be computationally burdensome and therefore k(100) or perhaps k(500) may be preferred if T or NT
is large.
ALGorithm(string): This option specifies the algorithm used for prediction. If not specified, simplex projection (a locally weighted average) is used. Valid options include simplex and smap, the latter of which is a sequential locally weighted global linear mapping (or S-map as noted previously). In the case of the xmap subcommand where two variables predict each other, the algorithm(smap) invokes something analogous to a distributed lag model with E + 1 predictors (including a constant term c) and, thus, E + 1 locally-weighted coefficients for each predicted observation/target vector—because each predicted observation has its own type of regression done with k neighbours as rows and E + 1 coefficients as columns. As noted below, in this case special options are available to save these coefficients for post-processing but, again, it is not actually a regression model and instead should be seen as a manifold.
RANDomize: When splitting the observations into library and prediction sets, by default the oldest observations go into the library set and the newest observations to the prediction set. Though if the randomize option is specified, the data is allocated into the two sets in a random fashion. If the replicate option is specified, then this randomization is enabled automatically.
REPlicate(integer): The explore subcommand uses a random 50/50 split for simplex projection and S-maps, whereas the xmap subcommand selects the observations randomly for library construction if the size of the library L is smaller than the size of all available observations. In these cases, results may be different in each run because the embedding vectors (i.e., the E-dimensional points) used to reconstruct a manifold are chosen at random. The replicate option takes advantages of this to allow repeating the randomization process and calculating results each time. This is akin to a nonparametric bootstrap without replacement, and is commonly used for inference using confidence intervals in EDM (Tsonis et al., 2015; van Nes et al., 2015; Ye et al., 2015b). When replicate is specified, such as replicate(50), mean values and the standard deviations of the results are reported across the 50 runs by default. As we note below, it is possible to save all estimates for post-processing using typical Stata commands such as svmat, allowing the graphing of results or finding percentile-based with the pctile command.
PREDICTionsave(variable): This option allows you to save the edm predictions as a variable, which could be useful for plotting and diagnosis.
COPREDICTionsave(variable): This option allows you to save the copredictions as a variable. You must specify the copredictvar(variables) options for this to work.
copredictvar(variable): This option specifies the variable used for coprediction. A second prediction is run for each configuration of E, library, etc., using the same library set but with a prediction set built from the lagged embedding of this variable.
Predictionhorizon(integer): This option adjusts the default number of observations ahead which we predict. By default, the explore mode predict τ observations ahead and the xmap mode uses p(0). This parameter can be negative.
DETails: By default, only mean values and standard deviations are reported when the replicate option is specified. The details option overrides this behaviour by providing results for each individual run. Irrespective of using this option, all results can be saved for post-processing.
CI(integer): When used with replicate() or crossfold(), this option reports the confidence interval for the mean of the estimates (MAE and/or ρ), as well as the percentiles of their distribution. The first row of output labelled “Est. mean CI” reports the estimated confidence interval of the mean ρ, assuming that ρ has a normal distribution—estimated as the corrected sample standard deviation (with N-1 in the denominator) divided by the squared root of the number of replications. The reported range can be used to compare mean ρ across different (hyper) parameter values (e.g., different E, θ, or L) using the same datasets as if the sample was the entire population (such that uncertainty is reduced to 0 when the number of replications →∞). These intervals can be used to test which (hyper) parameter values best describe a sample, as might be typically used when using crossfold validation methods. The row labelled with “Pc (Est.)” follows the same normality assumption and reports the estimated percentile values based on the corrected sample standard deviation of the replicated estimates. The row labelled “Pc (Obs.)” reports the actual observed percentile values from the replicated estimates. In both of these latter cases the percentile values offer alternative metrics for comparisons across distributions, which would be more useful for testing typical hypotheses about population differences in estimates across different (hyper) parameter values (e.g., different E, θ, or L), such as testing whether a dynamical system appears to be nonlinear in a population (i.e., testing whether ρ is maximized when θ > 0). The number specified within the ci() bracket determines the confidence level and the locations of the percentile cut-offs. For example, ci(90) instructs edm to return 90% CI as well as the cut-off values for the 5th and 95th percentile values—because ρ and MAE values cannot or are not expected to take on negative values, we typically prefer one-tailed hypothesis tests and therefore would use ci(90) to get a one-tailed 95% interval. These estimated ranges are also included in the e() return list as a series of scalars with names starting with “ub” for upper bound and “lb” for lower bound values of the CIs. These return values can be used for further post-processing.
seed(integer): This option specifies the seed used for the random number. In some
special cases users may wish to use this in order to keep library and prediction sets the same
across simplex projection and S-mapping with a single variable, or across multiple CCM runs with
different variables. Note: if set rng has been used to change Stata's random number generation
algorithm, then edm will temporarily change it the default 64 bit Mersenne twister internally.
strict: When this option is specified, the computation will fail if the requested number of neighboring observations is not present. By default, if not all of the request neighbors are available, the computation will just continue using as many as possible.
ALLOWMISSing This option allows observations with missing values to be used in the manifold. Vectors with at least one non-missing values will be used in the manifold construction. Distance computations are adapted to allow missing values when this option is specified.
MISSINGdistance(real): This option allows users to specify the assumed distance between missing values and any values (including missing) when estimating the Euclidean distance of the vector. This enables computations with missing values. The option implies allowmissing. By default, the distance is set to the expected distance of two random draws in a normal distribution, which equals to 2/sqrt(pi) * standard deviation of the mapping variable.
EXTRAembed(variables): This option allows incorporating additional variables into the
embedding (multivariate embedding), e.g. extra(z l.z). Time series lists are unabbreviated here,
e.g. extra(L(1/3).z) will be equivalent to extra(L1.z L2.z L3.z). Normally, lagged versions of the
extra variables are not included in the embedding, however the syntax extra(z(e)) includes E lags
of z in the embedding.
dt: This option allows automatic inclusion of the timestamp differencing in the
embedding. There will be E dt variables included for an embedding with E dimensions. By default,
the weights used for these additional variables equal to the standard deviation of the main
mapping variable divided by the standard deviation of the time difference. This can be overridden by
the dtweight() option. dt option will be ignored when running with data with no sampling
variation in the time lags. The first dt variable embeds the time of the between the most recent
observation and the time of the corresponding target/predictand.
reldt: This option, to be read as 'relative dt', is like the 'dt' option above in that it includes E extra variables for an embedding with E dimensions. However the timestamp differences added are not the time between the corresponding observations, but the time of the target/predictand minus the time of the lagged observations.
DTWeight(real): This option specifies the weight used for the timestamp differencing variable.
DTSave(variable): This option allows users to save the internally generated timestamp differencing variable.
nthreads(integer): The number of threads the C++ plugin will use for parallel computations. The default is the number of cores available on the host computer.
idw(real): This parameter is used when xtset indicates that the current dataset is panel data. Then, idw specifies a penalty that is added to the distances between points in the manifold which correspond to observations from different panels. By default idw is 0, so the data from all panels is mixed together and treatly equally. If idw(-1) is set (or any other negative value), then the weight is treated as 'infinity', so neighbours will never be selected which cross the boundaries between panels. Setting idw(-1) with k(-1) means we may use a different number of neighbors for different predictions (i.e. if the panels are unbalanced).
lowmemory: It is possible that RAM may be depleted while running edm on large datasets with large values of E. The lowmemory flag directs the plugin to try to save as much space as possible by more efficiently using memory, though for small datasets this will likely slow down the computations by a small but noticeable amount.
Besides the shared parameters, edm explore supports the following extra options:
CROSSfold(integer): This option asks the program to run a cross-fold validation of the
predicted variables. crossfold(5) indicates a 5-fold cross validation. Note that this cannot be used
together with replicate. This option is only available with the explore subcommand.
full: When this option is specified, the explore command will use all possible observations in the manifold construction instead of the default 50/50 split. This is effectively the same as leave-one-out cross-validation as the observation itself is not used for the prediction.
reportrawe: By default, the program reports the actual E used in the manifold. With this option, the program will only report the number of dimensions constructed from the main variable.
Library(numlist ascending): This option specifies the total library size L used for
the manifold reconstruction. Varying the library size is used to estimate the convergence property
of the cross-mapping, with a minimum value Lmin = E + 2 and the maximum equal to the total number of
observations minus sufficient lags (e.g., in the time-series case without missing data this is Lmax
= T + 1 – E). An error message is given if the L value is beyond the allowed range. To assess the
rate of convergence (i.e., the rate at which ρ increases as L grows), the full range of library
sizes at small values of L can be used, such as if E = 2 and T = 100, with the setting then perhaps
being library(4(1)25 30(5)50 54(15)99). This option is only available with the xmap subcommand.
SAVEsmap(string): This option allows smap coefficients to be stored in variables with
a specified prefix. For example, specifying edm xmap x y, algorithm(smap) savesmap(beta) k(-1) will
create a set of new variables such as beta1_b0_rep1. The string prefix (e.g., 'beta') must not be
shared with any variables in the dataset, and the option is only valid if the algorithm(smap) is
specified. In terms of the saved variables such as beta1_b0_rep1, the first number immediately after
the prefix ‘beta’ is 1 or 2 and indicates which of the two listed variables is treated as the
dependent variable in the cross-mapping (i.e., the direction of the mapping). For the edm xmap x y
case, variables starting with beta1_ contain coefficients derived from the manifold M_X created
using the lags of the first variable ‘x’ to predict Y, or Y|M_X. This set of variables therefore
store the coefficients related to ‘x’ as an outcome rather than a predictor in CCM. Keep in mind
that any Y→X effect associated with the beta1_ prefix is shown as Y|M_X, because the outcome is used
to cross-map the predictor, and thus the reported coefficients will be scaled in the opposite
direction of a typical regression (because in CCM the outcome variable predicts the cause). To get
more familiar regression coefficients (which will be locally weighted), variables starting with
beta2_ store the coefficients estimated in the other direction, where the second listed variable ‘y’
is used for the manifold reconstruction M_Y for the mapping X|M_Y in the “edm xmap x y” case, testing
the opposite X→Y effect in CCM, but with reported S-map coefficients that map to a Y→X regression.
We appreciate that this may be unintuitive, but because CCM causation is tested by predicting the
causal variable with the outcome, to get more familiar regression coefficients requires reversing
CCM’s causal direction to a more typical predictor→outcome regression logic. This can be clarified
by reverting to the conditional notation such as X|M_Y, which in CCM implies a left-to-right X→Y
effect, but for the S-map coefficients will be scaled as a locally-weighted regression in the
opposite direction Y→X. Moving on, following the 1 and 2 is the letter b and a number. The numerical
labeling scheme generally follows the order of the lag for the main variable and then the order of
the extra variables introduced in the case of multivariate embedding. b0 is a special case which
records the coefficient of the constant term in the regression. The final term rep1 indicates the
coefficients are from the first round of replication (if the replicate() option is not used then
there is only one). Finally, the coefficients are saved to match the observation t in the dataset
that is being predicted, which allows plotting each of the E estimated coefficients against time
and/or the values of the variable being predicted. The variables are also automatically labelled for
clarity. This option is only available with the xmap subcommand.
DIrection(string): This option allows users to control whether the cross mapping is
calculated bidirectionally or unidirectionally, the latter of which reduces computation times if
bidirectional mappings are not required. Valid options include “oneway” and “both”, the latter of
which is the default and computes both possible cross-mappings. When oneway is chosen, the first
variable listed after the xmap subcommand is treated as the potential dependent variable following
the conventions in the regression syntax of Stata such as the‘reg’ command, so edm xmap x y, direction(oneway) produces the cross-mapping Y|M_X, which pertains to a Y→X effect. This is
consistent with the beta1_ coefficients from the previous savesmap(beta) option. On this point, the
direction(oneway) option may be especially useful when an initial “edm xmap x y”procedure shows
convergence only for a cross-mapping Y|M_X, which pertains to a Y→X effect. To save time with large
datasets, any follow-up analyses with the algorithm(smap) option can then be conducted with edm xmap x y, algorithm(smap) savesmap(beta) direction(oneway). To make this easier there is also a
simplified oneway option that implies direction(oneway). This option is only available with the
xmap subcommand.
oneway: This option is equivalent to direction(oneway)
Options for update subcommand
The update subcommand supports the following options:
develop: This option updates the command to its latest development version. The development version usually contains more features but may be less tested compared with the older version distributed on SSC.
replace: This option specifies whether you allow the update to override your local ado files.
Examples
Chicago crime dataset example (included in the auxiliary file)
Updates
To install the stable version or upgrade directly through Stata:
To install the development version directly through Stata: