Stable law estimation functions

Description

A collection of methods to estimate the four parameters of stable laws. The package also provides functions to compute the characteristic function and tools to run Monte Carlo simulations.

Details

The main functions of the package are briefly described below:

main function:

Estim is the most useful function of the package. It estimates of the parameters and the asymptotic properties of the estimators.

estimation function:

the methods provided so far are the maximum-likelihood (MLParametersEstim), the generalised method of moment with finite (GMMParametersEstim) or continuum (CgmmParametersEstim) moment conditions, the iterative Koutrouvelis regression method (KoutParametersEstim) and the fast Kogon-McCulloch method used for first guess estimation (IGParametersEstim).

characteristic function:

the characteristic function (ComplexCF) and its Jacobian (jacobianComplexCF) can be computed and will return a vector (respectively a matrix) of complex numbers.

Monte Carlo simulation

Estim_Simulation is a tool to run Monte Carlo simulations with flexible options to select the estimation method, the Monte Carlo control parameters, compute statistical summaries or save results to a file.

Note

Version 1 of this package had a somewhat restricted license since it needed package akima in some computations.

In version 2 of the package we implemented a 2D interpolation routine and removed the dependency on akima. Therefore, StableEstim is now under GPL license. The package is related to upcoming work by the authors where the different methods are compared using MC simulations.

Author(s)

Tarak Kharrat, Georgi N. Boshnakov

References

Carrasco M and Florens J (2000). “Generalization of GMM to a continuum of moment conditions.” Econometric Theory, 16(06), pp. 797–834.

Carrasco M and Florens J (2002). “Efficient GMM estimation using the empirical characteristic function.” IDEI Working Paper, 140.

Carrasco M and Florens J (2003). “On the asymptotic efficiency of GMM.” IDEI Working Paper, 173.

Carrasco M, Chernov M, Florens J and Ghysels E (2007). “Efficient estimation of general dynamic models with a continuum of moment conditions.” Journal of Econometrics, 140(2), pp. 529–573.

Carrasco M, Florens J and Renault E (2007). “Linear inverse problems in structural econometrics estimation based on spectral decomposition and regularization.” Handbook of econometrics, 6, pp. 5633–5751.

Carrasco M and Kotchoni R (2010). “Efficient estimation using the characteristic function.” Mimeo. University of Montreal.

Nolan J (2001). “Maximum likelihood estimation and diagnostics for stable distributions.” L'evy processes: theory and applications, pp. 379–400.

Nolan JP (2012). Stable Distributions - Models for Heavy Tailed Data. Birkhauser, Boston. In progress, Chapter 1 online at academic2.american.edu/\simjpnolan.

Hansen LP (1982). “Large sample properties of generalized method of moments estimators.” Econometrica: Journal of the Econometric Society, pp. 1029–1054.

Hansen LP, Heaton J and Yaron A (1996). “Finite-sample properties of some alternative GMM estimators.” Journal of Business & Economic Statistics, 14(3), pp. 262–280.

Feuerverger A and McDunnough P (1981). “On efficient inference in symmetric stable laws and processes.” Statistics and Related Topics, 99, pp. 109–112.

Feuerverger A and McDunnough P (1981). “On some Fourier methods for inference.” Journal of the American Statistical Association, 76(374), pp. 379–387.

Schmidt P (1982). “An improved version of the Quandt-Ramsey MGF estimator for mixtures of normal distributions and switching regressions.” Econometrica: Journal of the Econometric Society, pp. 501–516.

Besbeas P and Morgan B (2008). “Improved estimation of the stable laws.” Statistics and Computing, 18(2), pp. 219–231.

See Also

fBasics:::.mleStableFit, fBasics:::.qStableFit

package stabledist

Class "Best_t"

Description

Class used to store the result of function ComputeBest_t.

Objects from the Class

Objects can be created by calls of the form new("Best_t", theta, nbt, tvec, detVal, convcode, ...), where the user can specify some/all of the inputs or call function ComputeBest_t.

Slots

theta:

Object of class "vector"; values of the 4 parameters.

nbt:

Object of class "vector"; number of points used in the minimisation.

tvec:

Object of class "list"; values of the best t-vectors.

detVal:

Object of class "vector"; values of the optimal determinant found after minimisation.

convcode:

Convergence code.

Methods

+

signature(e1 = "Best_t", e2 = "Best_t"): sum objects from class Best_t.

initialize

signature(.Object = "Best_t"): initialise an object from class Best_t as described above.

show

signature(object = "Best_t"): print a summary of the object.

See Also

ComputeBest_t

Estimate parameters of stable laws using a Cgmm method

Description

Estimate the four parameters of stable laws using generalised method of moments based on a continuum of complex moment conditions (Cgmm) due to Carrasco and Florens. Those moments are computed by matching the characteristic function with its sample counterpart. The resulting (ill-posed) estimation problem is solved by a regularisation technique.

Usage

CgmmParametersEstim(x, type = c("2S", "IT", "Cue"), alphaReg = 0.01,
                    subdivisions = 50,
                    IntegrationMethod = c("Uniform", "Simpson"),
                    randomIntegrationLaw = c("unif", "norm"),
                    s_min = 0, s_max = 1,
                    theta0 = NULL,
                    IterationControl = list(),
                    pm = 0, PrintTime = FALSE,...)

Arguments

Details

The moment conditions The moment conditions are given by:

g_t(X,\theta)=g(t,X;\theta)= e^{itX} - \phi_{\theta}(t)

If one has a sample x_1,\dots,x_n of i.i.d realisations of the same random variable X, then:

\hat{g}_n(t,\theta) = \frac{1}{n}\sum_{i=1}^n g(t,x_i;\theta) = \phi_n(t) -\phi_\theta(t),

where \phi_n(t) is the eCF associated with the sample x_1,\dots,x_n, defined by \phi_n(t)= \frac{1}{n} \sum_{j=1}^n e^{itX_j}. Objective function

Following Carrasco et al. (2007), Proposition 3.4, the objective function to minimise is given by:

obj(\theta)=\overline{\underline{v}^{\prime}}(\theta)[\alpha_{Reg} \mathcal{I}_n+C^2]^{-1}\underline{v}(\theta)

where:

\underline{v} = [v_1,\ldots,v_n]^{\prime};

v_i(\theta) = \int_I \overline{g_i}(t;\hat{\theta}^1_n) \hat{g}(t;\theta) \pi(t) dt.

I_n

is the identity matrix of size n.

C

is a n \times n matrix with (i,j)th element given by c_{ij} = \frac{1}{n-4}\int_I \overline{g_i}(t;\hat{\theta}^1_n) g_j(t;\hat{\theta}^1_n) \pi(t) dt.

To compute C and v_i() we will use the function IntegrateRandomVectorsProduct. The IterationControl If type = "IT" or type = "Cue", the user can control each iteration using argument IterationControl, which should be a list which contains the following elements:

NbIter:

maximum number of iterations. The loop stops when NBIter is reached; default = 10.

PrintIterlogical:

if set to TRUE the values of the current parameter estimates are printed to the screen at each iteration; default = TRUE.

RelativeErrMax:

the loop stops if the relative error between two consecutive estimation steps is smaller then RelativeErrMax; default = 1e-3.

Value

a list with the following elements:

Note

nlminb as used to minimise the Cgmm objective function.

References

Carrasco M, Florens J (2000). “Generalization of GMM to a continuum of moment conditions.” Econometric Theory, 16(06), 797–834.

Carrasco M, Florens J (2002). “Efficient GMM estimation using the empirical characteristic function.” IDEI Working Paper, 140.

Carrasco M, Florens J (2003). “On the asymptotic efficiency of GMM.” IDEI Working Paper, 173.

Carrasco M, Chernov M, Florens J, Ghysels E (2007). “Efficient estimation of general dynamic models with a continuum of moment conditions.” Journal of Econometrics, 140(2), 529–573.

Carrasco M, Kotchoni R (2010). “Efficient estimation using the characteristic function.” Mimeo. University of Montreal.

See Also

Estim, GMMParametersEstim, IntegrateRandomVectorsProduct

Examples

## general inputs
theta <- c(1.45, 0.55, 1, 0)
pm <- 0
set.seed(2345)
x <- rstable(50, theta[1], theta[2], theta[3], theta[4], pm)

## GMM specific params
alphaReg <- 0.01
subdivisions <- 20
randomIntegrationLaw <- "unif"
IntegrationMethod <- "Uniform"

## Estimation
twoS <- CgmmParametersEstim(x = x, type = "2S", alphaReg = alphaReg, 
                          subdivisions = subdivisions, 
                          IntegrationMethod = IntegrationMethod, 
                          randomIntegrationLaw = randomIntegrationLaw, 
                          s_min = 0, s_max = 1, theta0 = NULL, 
                          pm = pm, PrintTime = TRUE)
twoS

Compute the characteristic function of stable laws

Description

Theoretical characteristic function (CF) of stable laws under parametrisation ‘S0’ or ‘S1’. See Nolan (2013) for more details.

Usage

ComplexCF(t, theta, pm = 0)

Arguments

Details

For more details about the different parametrisation of the CF, see Nolan(2012).

Value

vector of complex numbers with dimension length(t).

References

Nolan JP (2012). Stable Distributions - Models for Heavy Tailed Data. Birkhauser, Boston. In progress, Chapter 1 online at academic2.american.edu/\simjpnolan.

See Also

jacobianComplexCF

Examples

## define the parameters
nt <- 10
t <- seq(0.1, 3, length.out = nt)
theta <- c(1.5, 0.5, 1, 0)
pm <- 0

## Compute the characteristic function
CF <- ComplexCF(t = t, theta = theta, pm = pm)
CF

Monte Carlo simulation to investigate the optimal number of points to use in the moment conditions

Description

Runs Monte Carlo simulation for different values of \alpha and \beta and computes a specified number of t-points that minimises the determinant of the asymptotic covariance matrix.

Usage

ComputeBest_t(AlphaBetaMatrix = abMat, nb_ts = seq(10, 100, 10),
              alphaReg = 0.001, FastOptim = TRUE, ...)

Arguments

Value

a list containing slots from class Best_t-class corresponding to one value of the parameters \alpha and \beta.

See Also

ComputeBest_tau, Best_t-class

Run Monte Carlo simulation to investigate the optimal \tau

Description

Runs Monte Carlo simulation to investigate the optimal number of points to use when one of the reduced spacing schemes is considered.

Usage

ComputeBest_tau(AlphaBetaMatrix = abMat, nb_ts = seq(10, 100, 10),
                tScheme = c("uniformOpt", "ArithOpt"),
                Constrained = TRUE, alphaReg = 0.001, ...)

Arguments

Value

a list containing slots from class Best_t-class corresponding to one value of the parameters \alpha and \beta.

See Also

ComputeBest_t, Best_t-class

Duration

Description

Compute the duration between 2 time points.

Usage

ComputeDuration(t_init, t_final, OneNumber = FALSE)

Arguments

Value

a numeric of length 1 or 3 depending on the value of OneNumber flag.

See Also

PrintDuration, PrintEstimatedRemainingTime.

Examples

ti <- getTime_()
for (i in 1:100) x <- i*22.1
tf <- getTime_()
ComputeDuration(ti,tf)

First root of the empirical characteristic function

Description

Computes the first root of the real part of the empirical characteristic function.

Usage

ComputeFirstRootRealeCF(x, ..., tol = 0.001, maxIter = 100,
                        lowerBand = 1e-04, upperBand = 30)

Arguments

Details

The Welsh algorithm is first applied. If it fails to provide a satisfactory value (< tol), a graphical/ numerical approach is used. We first plot the real part of the eCF vs t in order to determine the first zero directly and use it as the initial guess of a numerical minimisation routine.

Value

numeric: first zero of the real part of the eCF.

References

Welsh AH (1986). “Implementing empirical characteristic function procedures.” Statistics & probability letters, 4(2), 65–67.

See Also

ComplexCF

Examples

set.seed(345)
x <- rstable(500, 1.5, 0.5)
ComputeFirstRootRealeCF(x)

Parse an output file to create a summary object (list)

Description

Parses the file saved by Estim_Simulation and re-creates a summary list identical to the one produced by Estim_Simulation when StatSummary is set to TRUE.

Usage

ComputeStatObjectFromFiles(files, sep_ = ",",
                           FctsToApply = StatFcts,
                           headers_=TRUE,readSizeFrom=1,
                           CheckMat=TRUE,
                           tolFailCheck=tolFailure,
                           MCparam=1000,...)

Arguments

Details

The same sample sizes are assumed for all the files and we also assume a different set of parameters (alpha,beta) within each file (one and one only).

This function is particularly useful when simulations are run in parallel on different computers/CPUs and the output files are collected afterwards. This function is also used to create the Latex summary table: see TexSummary.

Some examples are provided in the example folder.

Value

a list of length 4 containing a summary matrix object associated to each parameter.

See Also

Estim_Simulation

Concatenates output files.

Description

Creates a unique file by concatenating several output files associated to one set of parameters.

Usage

ConcatFiles(files, sep_ = ",", outfile, headers_ = TRUE,
            DeleteIfExists=TRUE)

Arguments

Details

The files to be concatenated should be related to the same set of parameters alpha and beta. The function stops if one of the file contains 2 (or more) different set of parameters (the function compares the values of columns 1 and 2 row by row) or if the set of parameters within one file is different from the one from other files.

Value

Returns an output file outfile saved in the working directory.

See Also

Estim_Simulation

Estimate parameters of stable laws

Description

Estimates the four parameters of stable distributions using one of the methods implemented in StableEstim. This is the main user-level function but the individul methods are available also as separate functions.

Usage

Estim(EstimMethod = c("ML", "GMM", "Cgmm","Kout"), data, theta0 = NULL,
      ComputeCov = FALSE, HandleError = TRUE, ...)

Arguments

Details

Estim is the main estimation function in package StableEstim.

This function should be used in priority for estimation purpose as it provides more information about the estimator. However, user needs to pass the appropriate parameters to the selected method in .... See the documentation of the selected method.

Asymptotic Confidence Intervals: The normal asymptotic confidence intervals (CI) are computed. The user can set the level of confidence by inputing the level argument (in the "\dots"); default level=0.95. The theoretical justification for asymptotic normal CI can be found in the references for the individual methods. Note the CI's are not computed for the Koutrouvelis regression method.

Value

an object of class Estim, see Estim-class for more details

See Also

CgmmParametersEstim, GMMParametersEstim, MLParametersEstim, KoutParametersEstim for the individual estimation methods;

IGParametersEstim for fast computation of initial values.

Examples

## general inputs
theta <- c(1.45, 0.55, 1, 0)
pm <- 0
set.seed(2345)
x <- rstable(200, theta[1], theta[2], theta[3], theta[4], pm)

objKout <- Estim(EstimMethod = "Kout", data = x, pm = pm, 
                     ComputeCov = FALSE, HandleError = FALSE, 
                     spacing = "Kout")

Class "Estim"

Description

Class for storing the results of estimating parameters of stable laws, output of function Estim().

Objects from the Class

Objects can be created by calls of the form new("Estim", par, ...). Users can provide some (or all) of the inputs stated below to create an object from this class or call function Estim with appropriate arguments.

Slots

par:

numeric(4), values of the 4 estimated parameters.

par0:

numeric(4), initial values for the 4 parameters.

vcov:

object of class "matrix" (4 x 4), representing the covariance matrix of the estimated parameters.

confint:

object of class "matrix" (4 x 4), representing the confidence interval computed at a specific level (attribute of the object).

data:

numeric(), the data used to compute the estimates.

sampleSize:

numeric(1), length of the data.

others:

list(), further information about the estimation method.

duration:

numeric(1), duration in seconds.

failure:

numeric(1), represents the status of the procedure: 0 failure or 1 success.

method:

Object of class "character", description of the parameters used in the estimation.

Methods

initialize

signature(.Object = "Estim"): creates an object of this class using the inputs described above.

show

signature(object = "Estim"): summarised print of the object.

See Also

Estim

Monte Carlo simulation

Description

Runs Monte Carlo simulation for a selected estimation method. The function can save a file and produce a statistical summary.

Usage

Estim_Simulation(AlphaBetaMatrix = abMat, SampleSizes = c(200, 1600),
                 MCparam = 100, Estimfct = c("ML", "GMM", "Cgmm","Kout"),
                 HandleError = TRUE, FctsToApply = StatFcts,
                 saveOutput = TRUE, StatSummary = FALSE,
                 CheckMat = TRUE, tolFailCheck = tolFailure,
                 SeedOptions=NULL, ...)

Arguments

Details

Error Handling It is advisable to set it to TRUE when the user is planning to launch long simulations as it will prevent the procedure from stopping if an error occurs for one sample data. The estimation function will produce a vector of NA as estimated parameters related to this (error generating) sample data and move on to the next Monte Carlo step. Statistical summary The function is able to produce a statistical summary of the Monte Carlo simulation for each parameter (slices of the list). Each slice is a matrix where the rows represents the true values of the parameters and the columns the statistical information. In all cases, the following quantities are computed:

sample size:

the sample size used to produce the simulated data.

alphaT, betaT:

the true values of the parameters.

failure:

the number of times the procedure failed to produce relevant estimation.

time:

the average running time in seconds of the estimation procedure

Besides, the (vector of character) FctsToApply controls the other quantities to be computed by providing the name of the function object to be applied to the vector of estimated parameters. The signature of the function should be of the form fctName = function(p,...){...}, where p is the vector (length(p) = MCparam) of parameter estimates and ... is the extra arguments to be passed the function.

By default, the functions from StatFcts will be applied but the user can pass his own functions by providing their names in argument FctsToApply and their definitions in the global environment.

Note that if CheckMat is set to TRUE, the estimation is considered failed if the squared error (of the first 2 parameters alpha and beta) is larger than tolFailCheck.

Output file

Setting saveOutput to TRUE will have the side effect of saving a csv file in the working directory. This file will have MCparam * length(SampleSizes) lines and its columns will be:

alphaT, betaT:

the true values of the parameters.

data size:

the sample size used to generate the simulated data.

seed:

the seed value used to generate the simulated data.

alphaE, betaE, gammaE, deltaE:

the estimates of the 4 parameters.

failure:

binary: 0 for success, 1 for failure.

time:

estimation running time in seconds.

The file name is informative to let the user identify the values of the true parameters, the MC parameters, as well as the options selected for the estimation method.

The csv file is updated after each MC estimation, which is useful when the simulation stops before it finishes. Besides, using the check-pointing mechanism explained below, the simulation can re-start from where it stopped.

Check-pointing. Checkpointing is the act of saving enough program state and results so far calculated that a computation can be stopped and restarted. The way we did it here is to save a text file with some useful information about the state of the estimation. This text file is updated after each MC iteration and read at the beginning of function Estim_Simulation to allow the simulation to re-start from where it stopped. This file is deleted at the end of the simulation procedure.

SeedOptions. Users who do not want to control the seed generation can ignore this argument (its default value is NULL). This argument can be more useful when one wants to cut the simulation (even for one parameter value) into pieces. In that case, the user can control which part of the seed vector to use.

MCtot:

total values of MC simulations in the entire process.

seedStart:

starting index in the seed vector. The vector extracted will be of size MCparam.

Value

If StatSummary is set to TRUE, a list with 4 components (corresponding to the 4 parameters) is returned. Each component is a matrix. If SaveOutput is set to TRUE, only a csv file is saved and nothing is returned (if StatSummary is FALSE). If both are FALSE, the function stops.

See Also

Estim, CgmmParametersEstim, GMMParametersEstim, MLParametersEstim

Estimate parameters of stable laws using a GMM method

Description

Estimate parameters of stable laws using generalised method of moments (GMM) with finite number of moment conditions. It uses a regularisation technique to make the method more robust (when the number of moment condition is large) and allows different schemes to select where the moment conditions are computed.

Usage

GMMParametersEstim(x, algo = c("2SGMM", "ITGMM", "CueGMM"),
                   alphaReg = 0.01,
                   regularization = c("Tikhonov", "LF", "cut-off"),
                   WeightingMatrix = c("OptAsym", "DataVar", "Id"),
                   t_scheme = c("equally", "NonOptAr", "uniformOpt",
                                "ArithOpt", "VarOpt", "free"),
                   theta0 = NULL,
                   IterationControl = list(),
                   pm = 0, PrintTime = FALSE, ...)

Arguments

Details

The moment conditions

The moment conditions are given by:

g_t(X,\theta) = g(t,X;\theta)= e^{itX} - \phi_{\theta}(t)

If one has a sample x_1,\dots,x_n of i.i.d realisations of the same random variable X, then:

\hat{g}_n(t,\theta) = \frac{1}{n}\sum_{i=1}^n g(t,x_i;\theta) = \phi_n(t) -\phi_\theta(t),

where \phi_n(t) is the eCF associated to the sample x_1,\dots,x_n, and defined by \phi_n(t)= \frac{1}{n} \sum_{j=1}^n e^{itX_j}.

Objective function

obj{\theta} = < K^{-1/2} \hat{g}_n(.;\theta),K^{-1/2}\hat{g}_n(.;\theta)>,

where K^{-1}f denotes the solution \varphi (when it exists) of the equation K \varphi=f and K^{-1/2}=(K^{-1})^{1/2}. The optimal choice of the Weighting operator K (a matrix in the GMM case) and its estimation are discussed in Hansen (1982).

Weighting operator (Matrix)

OptAsym:

the optimal asymptotic choice as described by Hansen. The expression of the components of this matrix could be found for example in Feuerverger and McDunnough (1981b).

DataVar:

the covariance matrix of the data provided.

Id:

the identity matrix.

the t-scheme

One of the most important features of this method is that it allows the user to choose how to place the points where the moment conditions are evaluated. The general rule is that users can provide their own set of points (option "free") or choose one of the other schemes. In the latter case they need to specify the number of points nb_t in argument "\dots" and eventually the lower and upper limit (by setting Constrained to FALSE and providing min_t and max_t) in the non-optimised case. If one of the optimised cases is selected, setting Constrained to FALSE will not constrain the choice of \tau, see below. We mean by optimised set of point, the set that minimises the (determinant) of the asymptotic covariance matrix as suggested by Schmidt (1982) and Besbeas and Morgan (2008).

6 options have been implemented:

"equally":

equally placed points in [min_t,max_t]. When provided, user's min_t and max_t will be used (when Coinstrained = FALSE). Otherwise, eps and An will be used instead (where An is the first zero of the eCF).

"NonOptAr":

non optimal arithmetic placement: t_j = \frac{j(j+1)}{nbt(nbt+1)}(max-eps); j=1,\dots,nbt, where max is the upper band of the set of points selected as discussed before.

"uniformOpt":

uniform optimal placement: t_j=j \tau, j=1,\dots, nbt

"ArithOpt":

arithmetic optimal placement: t_j=j(j+1) \tau, j=1,\dots nbt

"Var Opt":

optimal variance placement as explained above.

"free":

user needs to pass his own set of points in "\dots".

For the "ArithOpt" and "uniformOpt" schemes, the function to minimise is seen as a function of the real parameter \tau instead of doing a vectorial optimisition as in the "Var Opt" case. In the latter case, one can choose between a fast (but less accurate) optimisation routine or a slow (but more accurate) one by setting the FastOptim flag to the desired value.

The IterationControl

If type = "IT" or type = "Cue" the user can control each iteration by setting up the list IterationControl which contains the following elements:

NbIter:

maximum number of iteration. The loop stops when NBIter is reached; default = 10.

PrintIterlogical:

if set to TRUE, the value of the current parameter estimation is printed to the screen at each iteration; default = TRUE.

RelativeErrMax:

the loop stops if the relative error between two consecutive estimation steps is smaller than RelativeErrMax; default = 1e-3.

Value

a list with the following elements:

Note

nlminb was used for the minimisation of the GMM objective funcion and to compute tau in the "uniformOpt" and "ArithOpt" schemes. In the "Var Opt" scheme, optim was preferred. All those routines have been selected after running different tests using the summary table produced by package optimx for comparing the performance of different optimisation methods.

References

Hansen LP (1982). “Large sample properties of generalized method of moments estimators.” Econometrica: Journal of the Econometric Society, pp. 1029–1054.

Hansen LP, Heaton J and Yaron A (1996). “Finite-sample properties of some alternative GMM estimators.” Journal of Business & Economic Statistics, 14(3), pp. 262–280.

Feuerverger A and McDunnough P (1981). “On efficient inference in symmetric stable laws and processes.” Statistics and Related Topics, 99, pp. 109–112.

Feuerverger A and McDunnough P (1981). “On some Fourier methods for inference.” Journal of the American Statistical Association, 76(374), pp. 379–387.

Schmidt P (1982). “An improved version of the Quandt-Ramsey MGF estimator for mixtures of normal distributions and switching regressions.” Econometrica: Journal of the Econometric Society, pp. 501–516.

Besbeas P and Morgan B (2008). “Improved estimation of the stable laws.” Statistics and Computing, 18(2), pp. 219–231.

See Also

Estim, CgmmParametersEstim

Examples

## General data
theta <- c(1.5, 0.5, 1, 0)
pm <- 0
set.seed(345);
x <- rstable(100, theta[1], theta[2], theta[3], theta[4], pm)
##---------------- 2S free ----------------
## method specific arguments
regularization <- "cut-off"
WeightingMatrix <- "OptAsym"
alphaReg <- 0.005
t_seq <- seq(0.1, 2, length.out = 12)

## If you are just interested by the value
## of the 4 estimated parameters
t_scheme = "free"
algo = "2SGMM"
     
suppressWarnings(GMMParametersEstim(
    x = x, algo = algo, alphaReg = alphaReg, 
    regularization = regularization, 
    WeightingMatrix = WeightingMatrix, 
    t_scheme = t_scheme, 
    pm = pm, PrintTime = TRUE, t_free = t_seq))

Estimate parameters of stable laws by Kogon and McCulloch methods

Description

Kogon regression method is used together with the McCulloch quantile method to provide initial estimates of parameters of stable distributions.

Usage

IGParametersEstim(x, pm = 0, ...)

Arguments

Details

The parameters \gamma and \delta are estimated using the McCulloch(1986) quantile method from fBasics. The data is rescaled using those estimates and used to perform the Kogon regression method to estimate \alpha and \beta.

Value

a vector of length 4 containing the estimates of the 4 parameters.

References

Kogon SM and Williams DB (1998). “Characteristic function based estimation of stable distribution parameters.” A practical guide to heavy tailed data, pp. 311–335. McCulloch JH (1986). “Simple consistent estimators of stable distribution parameters.” Communications in Statistics-Simulation and Computation, 15(4), pp. 1109–1136.

See Also

Estim, McCullochParametersEstim

Examples

x <- rstable(200, 1.2, 0.5, 1, 0, pm = 0)
IGParametersEstim(x, pm = 0)

Integral outer product of random vectors

Description

Computes the integral outer product of two possibly complex random vectors.

Usage

IntegrateRandomVectorsProduct(f_fct, X, g_fct, Y, s_min, s_max,
                              subdivisions = 50,
                              method = c("Uniform", "Simpson"),
                              randomIntegrationLaw = c("norm","unif"),
                              ...)

Arguments

Details

The function computes the nx \times ny matrix C = \int_{s_{min}}^{s_{max}} f_s(X) g_s(Y) \pi(s) ds, such as the one used in the objective function of the Cgmm method. This is essentially an outer product with with multiplication replaced by integration.

There is no function in R to compute vectorial integration and computing C element by element using integrate may be very slow when length(X) (or length(y)) is large.

The function allows complex vectors as its integrands.

Value

an nx \times ny matrix C with elements:

c_{ij} = \int_{s_{min}}^{s_{max}} f_s(X_i) g_s(Y_j) \pi(s) ds .

Examples

## Define the integrand
f_fct <- function(s, x) {
    sapply(X = x, FUN = sampleComplexCFMoment, t = s, theta = theta)
}
f_bar_fct <- function(s, x) Conj(f_fct(s, x))

## Function specific arguments
theta <- c(1.5, 0.5, 1, 0)
set.seed(345)
X <- rstable(3, 1.5, 0.5, 1, 0)
s_min <- 0;
s_max <- 2
numberIntegrationPoints <- 10
randomIntegrationLaw <- "norm"

Estim_Simpson <-
    IntegrateRandomVectorsProduct(f_fct, X, f_bar_fct, X, s_min, s_max, 
                                  numberIntegrationPoints, 
                                  "Simpson", randomIntegrationLaw)
Estim_Simpson

Iterative Koutrouvelis regression method

Description

Iterative Koutrouvelis regression method with different spacing schemes (points where the eCF is computed).

Usage

KoutParametersEstim(x, theta0 = NULL,
                    spacing = c("Kout", "UniformSpac", "ArithSpac", "free"),
                    pm = 0, tol = 0.05, NbIter = 10, PrintTime = FALSE, ...)

Arguments

Details

spacing

4 options for the spacing scheme are implemented as described above. In particular:

UniformSpac, ArithSpac:

The user can specify the number of points to choose in both regression by inputting nb_t and nb_u. Otherwise the Koutrouvelis table will be used to compte them.

free:

The user is expected to provide t_points and u_points otherwise the Kout scheme will be used.

Value

a list with the following elements:

References

Koutrouvelis IA (1980). “Regression-type estimation of the parameters of stable laws.” Journal of the American Statistical Association, 75(372), pp. 918–928.

Koutrouvelis IA (1981). “An iterative procedure for the estimation of the parameters of stable laws: An iterative procedure for the estimation.” Communications in Statistics-Simulation and Computation, 10(1), pp. 17–28.

See Also

Estim

Examples

pm <- 0
theta <- c(1.45, 0.5, 1.1, 0.4)
set.seed(1235)
x <- rstable(200, theta[1], theta[2], theta[3], theta[4], pm = pm)
theta0 <- theta - 0.1
spacing <- "Kout"

KoutParametersEstim(x = x, theta0 = theta0, 
                        spacing = spacing, pm = pm)

Maximum likelihood (ML) method

Description

Uses the numerical ML approach described by Nolan to estimate the 4 parameters of stable law. The method may be slow for large sample size due to the use of numerical optimisation routine.

Usage

MLParametersEstim(x, theta0 = NULL, pm = 0, PrintTime = FALSE, ...)

Arguments

Details

The function performs the minimisation of the numerical (-)log-density of stable laws computed by function dstable from package stabledist.

After testing several optimisation routines, we have found out that the "L-BFGS-B" algorithm performs better with the ML method (faster, more accurate).

Value

a list with the following elements:

References

Nolan J (2001). “Maximum likelihood estimation and diagnostics for stable distributions.” L'evy processes: theory and applications, pp. 379–400.

See Also

Estim

Examples

theta <- c(1.5, 0.4, 1, 0)
pm <- 0
## 50 points does not give accurate estimation
## but it makes estimation fast for installation purposes
## use at least 200 points to get decent results.
set.seed(1333)
x <- rstable(50, theta[1], theta[2], theta[3], theta[4], pm)

## This example takes > 30 sec hence commented out
## Not run: 
  ML <- MLParametersEstim(x = x, pm = pm, PrintTime = TRUE)

## End(Not run)
## see the Examples folder for more examples.

Quantile-based method

Description

McCulloch quantile-based method.

Usage

McCullochParametersEstim(x)

Arguments

Details

The code is a modified version of function .qStableFit from package fBasics.

Value

numeric of length 4, represening the value of the 4 parameters.

References

McCulloch JH (1986). “Simple consistent estimators of stable distribution parameters.” Communications in Statistics-Simulation and Computation, 15(4), pp. 1109–1136.

See Also

Estim, IGParametersEstim

Examples

set.seed(333)
x <- rstable(500, 1.3, 0.4, 1, 0)
McCullochParametersEstim(x)

Print duration

Description

Print duration in human readable format.

Usage

PrintDuration(t, CallingFct = "")

Arguments

Details

The duration will be printed in the format: hours/minutes/seconds.

Value

Prints a character to the screen.

Examples

ti <- getTime_()
for (i in 1:100) x = i*22.1
tf <- getTime_()
duration <- ComputeDuration(ti, tf)
PrintDuration(duration, "test")

Estimated remaining time

Description

Prints the estimated remaining time in a loop. Useful in Monte Carlo simulations.

Usage

PrintEstimatedRemainingTime(ActualIter, ActualIterStartTime, TotalIterNbr)

Arguments

Details

Called at the end of each Monte Carlo step, this function will compute the duration of the actual step, an estimate of the remaining MC loops duration and prints the result to the screen in a human readable format using function PrintDuration.

See Also

PrintDuration, ComputeDuration.

Regularised Inverse

Description

Regularised solution of the (ill-posed) problem K\phi = r where K is a n \times n matrix, r is a given vector of length n. Users can choose one of the 3 schemes described in Carrasco and Florens (2007).

Usage

RegularisedSol(Kn, alphaReg, r,
               regularization = c("Tikhonov", "LF", "cut-off"),
               ...)

Arguments

Details

Following Carrasco and Florens(2007), the regularised solution of the problem K \phi=r is given by :

\varphi_{\alpha_{reg}} = \sum_{j=1}^{n} q(\alpha_{reg},\mu_j)\frac{<r,\psi_j >}{\mu_j} \phi_j ,

where q is a (positive) real function with some regularity conditions and \mu,\phi,\psi the singular decomposition of the matrix K.

The regularization parameter defines the form of the function q. For example, the "Tikhonov" scheme defines q(\alpha_{reg},\mu) = \frac{\mu^2}{\alpha_{reg}+\mu^2}.

When the matrix K is symmetric, the singular decomposition is replaced by a spectral decomposition.

Value

the regularised solution, a vector of length n.

References

Carrasco M, Florens J and Renault E (2007). “Linear inverse problems in structural econometrics estimation based on spectral decomposition and regularization.” Handbook of econometrics, 6, pp. 5633–5751.

See Also

solve

Examples

## Adapted from R examples for Solve 
## We compare the result of the regularized sol to the expected solution

hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, "+")}

K_h8 <- hilbert(8);
r8 <- 1:8

alphaReg_robust <- 1e-4
Sa8_robust <- RegularisedSol(K_h8,alphaReg_robust,r8,"LF")

alphaReg_accurate <- 1e-10
Sa8_accurate <- RegularisedSol(K_h8,alphaReg_accurate,r8,"LF")

## when pre multiplied by K_h8, the expected solution is 1:8
## User can check the influence of the choice of alphaReg

Objects exported from other packages

Description

These objects are imported from other packages. Follow the links below to see their documentation.

rstable

rstable

Default functions used to produce the statistical summary

Description

Default functions used to produce the statistical summary in the Monte Carlo simulations.

Details

The functions are:

Mean

.mean <- function(p,...) mean(p)

Min

.min <- function(p,...) min(p)

Max

.max <- function(p,...) max(p)

Sn

.Sn <- function(p,n,...) sqrt(n)*sd(p)

MSE

.MSE <- function(p,paramT,...) (1/length(p))*sum((p-paramT)^2)

Std error

.st.err <- function(p,...) sd(p)/sqrt(length(p))

To change the statistical summary, provide functions with similar signatures and pass a character vector containing the function names to Estim_Simulation.

LaTeX summary

Description

Creates a TeX table from a summary object or a vector of files.

Usage

TexSummary(obj, files = NULL, sep_ = ",", FctsToApply = StatFcts,
           caption = "Statistical Summary", label = "Simtab",
           digits = 3, par_index = 1, MCparam = 1000, ...)

Arguments

Details

Accepted values for par_index are c(1,2,3,4) or c("alpha","beta","gamma","delta") or mixed.

Some examples are provided in the example folder.

Value

a list of length length(par_index) whose elements are objects from class Latex (produced by toLatex)

See Also

Estim_Simulation, ComputeStatObjectFromFiles, xtable

Test approximate equality

Description

Tests the approximate equality of 2 objects. Useful for running tests.

Usage

expect_almost_equal(x, y, tolExpect = 0.001)

Arguments

Details

This function works with the expect_that function from package testhat to test equality between 2 objects with a given tolerance. It is used particularly for testing functions output. See the CF examples in the Examples folder.

See Also

expect_that,testthat

Examples

x <- 1.1
y <- 1.5
expect_almost_equal(x, y, 1)      # passes
## expect_almost_equal(x, y, 0.3) # fails

Default functions used to produce the statistical summary

Description

Default functions used to produce the statistical summary in the Monte Carlo simulations.

Usage

get.StatFcts()

Value

The functions computed are:

Mean

.mean <- function(p,...) mean(p)

Min

.min <- function(p,...) min(p)

Max

.max <- function(p,...) max(p)

Sn

.Sn <- function(p,n,...) sqrt(n)*sd(p)

MSE

.MSE <- function(p,paramT,...) (1/length(p))*sum((p-paramT)^2)

Std error

.st.err <- function(p,...) sd(p)/sqrt(length(p))

Users can define their own summaries by defining functions with similar signatures and passing a character vector containing the functions' names to Estim_Simulation.

Default set of parameters to pass to Estim_Simulation

Description

Default set of parameters to pass to Estim_Simulation, inspired by the one used by Koutrevelis (1980) in his simulation procedure.

Usage

get.abMat()

Value

a 2-columns matrix containing a wide range of parameters \alpha and \beta covering the entire parameters space.

Read time

Description

Reads the time when the function is called.

Usage

getTime_()

Value

a numeric.

See Also

PrintDuration, PrintEstimatedRemainingTime, ComputeDuration

Examples

ti <- getTime_()

Jacobian of the characteristic function of stable laws

Description

Numeric jacobian of the characteristic function (CF) as a function of the parameter \theta evaluated at a specific (vector) point t and a given value \theta.

Usage

jacobianComplexCF(t, theta, pm = 0)

Arguments

Details

The numerical derivation is obtained by a call to the function jacobian from package numDeriv. We have set up its arguments by default and the user is not given the option to modify them.

Value

a matrix length(t) \times 4 of complex numbers.

See Also

ComplexCF

Examples

## define the parameters
nt <- 10
t <- seq(0.1, 3, length.out = nt)
theta <- c(1.5, 0.5, 1, 0)
pm <- 0

## Compute the jacobian of the characteristic function
jack_CF <- jacobianComplexCF(t = t, theta = theta, pm = pm)

Complex moment condition based on the characteristic function

Description

Computes the moment condition based on the characteristic function as a complex vector.

Usage

sampleComplexCFMoment(x, t, theta, pm = 0)

Arguments

Details

The moment conditions

The moment conditions are given by:

g_t(X,\theta) = g(t,X;\theta)= e^{itX} - \phi_{\theta}(t)

If one has a sample x_1,\dots,x_n of i.i.d realisations of the same random variable X, then:

\hat{g}_n(t,\theta) = \frac{1}{n}\sum_{i=1}^n g(t,x_i;\theta) = \phi_n(t) - \phi_\theta(t) ,

where \phi_n(t) is the eCF associated to the sample x_1,\dots,x_n, and defined by \phi_n(t) = \frac{1}{n} \sum_{j=1}^n e^{itX_j}.

The function compute the vector of difference between the eCF and the CF at a set of given point t.

Value

a complex vector of length(t).

See Also

ComplexCF, sampleRealCFMoment

Examples

## define the parameters
nt <- 10   
t <- seq(0.1, 3, length.out = nt)
theta <- c(1.5, 0.5, 1, 0)
pm <- 0

set.seed(222)
x <- rstable(200, theta[1], theta[2], theta[3], theta[4], pm)

## Compute the characteristic function
CFMC <- sampleComplexCFMoment(x = x, t = t, theta = theta, pm = pm)
CFMC

Real moment condition based on the characteristic function

Description

Computes the moment condition based on the characteristic function as a real vector.

Usage

sampleRealCFMoment(x, t, theta, pm = 0)

Arguments

Details

The moment conditions

The moment conditions are given by:

g_t(X,\theta) = g(t,X;\theta) = e^{itX} - \phi_{\theta}(t) .

If one has a sample x_1,\dots,x_n of i.i.d realisations of the same random variable X, then:

\hat{g}_n(t,\theta) = \frac{1}{n}\sum_{i=1}^n g(t,x_i;\theta) = \phi_n(t) -\phi_\theta(t) ,

where \phi_n(t) is the eCF associated with the sample x_1,\dots,x_n, and defined by \phi_n(t) = \frac{1}{n} \sum_{j=1}^n e^{itX_j}.

The function compute the vector of difference between the eCF and the CF at a set of given point t. If length(t) = n, the resulting vector will be of length = 2n, where the first n components will be the real part and the remaining the imaginary part.

Value

a vector of length 2 * length(t).

See Also

ComplexCF, sampleComplexCFMoment

Examples

## define the parameters
nt <- 10   
t <- seq(0.1, 3, length.out = nt)
theta <- c(1.5, 0.5, 1, 0)
pm <- 0

set.seed(222)
x <- rstable(200, theta[1], theta[2], theta[3], theta[4], pm)

# Compute the characteristic function
CFMR <- sampleRealCFMoment(x = x, t = t, theta = theta, pm = pm)
CFMR