, Stepan Mazur [aut]
, Mathias Ljungdahl
[ctb]| Title: | Simulations and Statistical Inference for Linear Fractional Stable Motions |
|---|---|
| Description: | Contains functions for simulating the linear fractional stable motion according to the algorithm developed by Mazur and Otryakhin <doi:10.32614/RJ-2020-008> based on the method from Stoev and Taqqu (2004) <doi:10.1142/S0218348X04002379>, as well as functions for estimation of parameters of these processes introduced by Mazur, Otryakhin and Podolskij (2018) <arXiv:1802.06373>, and also different related quantities. |
| Authors: | Dmitry Otryakhin [aut, cre]
, Stepan Mazur [aut]
, Mathias Ljungdahl
[ctb] |
| Maintainer: | Dmitry Otryakhin <[email protected]> |
| License: | GPL-3 |
| Version: | 1.1.2 |
| Built: | 2025-01-06 02:46:22 UTC |
| Source: | https://gitlab.com/dmitry_otryakhin/tools-for-parameter-estimation-of-the-linear-fractional-stable-motion |
Computes the corresponding function value from Mazur et al. 2018.
a_p(p)a_p(p)
p |
power, real number from (-1,1) |
Mazur S, Otryakhin D, Podolskij M (2020). “Estimation of the linear fractional stable motion.” Bernoulli, 26(1), 226–252. https://doi.org/10.3150/19-BEJ1124.
a_tilda triggers a_tilda_cpp which is written in C++ and essentially performs the computation of the value.
a_tilda(N, m, M, alpha, H)a_tilda(N, m, M, alpha, H)
N |
a number of points of the lfsm. |
m |
discretization. A number of points between two nearby motion points |
M |
truncation parameter. A number of points at which the integral representing the definition of lfsm is calculated. So, after M points back we consider the rest of the integral to be 0. |
alpha |
self-similarity parameter of alpha stable random motion. |
H |
Hurst parameter |
Stoev S, Taqqu MS (2004). “Simulation methods for linear fractional stable motion and FARIMA using the Fast Fourier Transform.” Fractals, 95(1), 95-121. https://doi.org/10.1142/S0218348X04002379.
Defined for the two frequencies as
alpha_hat(t1, t2, k, path, H, freq)alpha_hat(t1, t2, k, path, H, freq)
t1, t2
|
real number such that t2 > t1 > 0 |
k |
increment order |
path |
sample path of lfsm on which the inference is to be performed |
H |
Hurst parameter |
freq |
Frequency of the motion. It can take two values: "H" for high frequency and "L" for the low frequency setting. |
The function triggers function phi, thus Hurst parameter is required only in high frequency case. In the low frequency, there is no need to assign H a value because it will not be evaluated.
Mazur S, Otryakhin D, Podolskij M (2020). “Estimation of the linear fractional stable motion.” Bernoulli, 26(1), 226–252. https://doi.org/10.3150/19-BEJ1124.
m<-45; M<-60; N<-2^14-M alpha<-1.8; H<-0.8; sigma<-0.3 freq='H' r=1; k=2; p=0.4; t1=1; t2=2 # Estimating alpha in the high frequency case # using preliminary estimation of H lfsm<-path(N=N,m=m,M=M,alpha=alpha,H=H, sigma=sigma,freq='L',disable_X=FALSE,seed=3)$lfsm H_est<-H_hat(p=p,k=k,path=lfsm) H_est alpha_est<-alpha_hat(t1=t1,t2=t2,k=k,path=lfsm,H=H_est,freq=freq) alpha_estm<-45; M<-60; N<-2^14-M alpha<-1.8; H<-0.8; sigma<-0.3 freq='H' r=1; k=2; p=0.4; t1=1; t2=2 # Estimating alpha in the high frequency case # using preliminary estimation of H lfsm<-path(N=N,m=m,M=M,alpha=alpha,H=H, sigma=sigma,freq='L',disable_X=FALSE,seed=3)$lfsm H_est<-H_hat(p=p,k=k,path=lfsm) H_est alpha_est<-alpha_hat(t1=t1,t2=t2,k=k,path=lfsm,H=H_est,freq=freq) alpha_est
Parameter freq is preserved to allow for investigation of the inference procedure in high frequency case.
ContinEstim(t1, t2, p, k, path, freq)ContinEstim(t1, t2, p, k, path, freq)
t1, t2
|
real number such that t2 > t1 > 0 |
p |
power |
k |
increment order |
path |
sample path of lfsm on which the inference is to be performed |
freq |
Frequency of the motion. It can take two values: "H" for high frequency and "L" for the low frequency setting. |
Mazur S, Otryakhin D, Podolskij M (2020). “Estimation of the linear fractional stable motion.” Bernoulli, 26(1), 226–252. https://doi.org/10.3150/19-BEJ1124.
m<-45; M<-60; N<-2^10-M alpha<-0.8; H<-0.8; sigma<-0.3 p<-0.3; k=3; t1=1; t2=2 lfsm<-path(N=N,m=m,M=M,alpha=alpha,H=H, sigma=sigma,freq='L',disable_X=FALSE,seed=3)$lfsm ContinEstim(t1,t2,p,k,path=lfsm,freq='L')m<-45; M<-60; N<-2^10-M alpha<-0.8; H<-0.8; sigma<-0.3 p<-0.3; k=3; t1=1; t2=2 lfsm<-path(N=N,m=m,M=M,alpha=alpha,H=H, sigma=sigma,freq='L',disable_X=FALSE,seed=3)$lfsm ContinEstim(t1,t2,p,k,path=lfsm,freq='L')
General estimation procedure for high frequency case when 1/alpha is not a natural number. "Unnecessary" parameter freq is preserved to allow for investigation of the inference procedure in low frequency case
GenHighEstim(p, p_prime, path, freq, low_bound = 0.01, up_bound = 4)GenHighEstim(p, p_prime, path, freq, low_bound = 0.01, up_bound = 4)
p |
power |
p_prime |
power |
path |
sample path of lfsm on which the inference is to be performed |
freq |
Frequency of the motion. It can take two values: "H" for high frequency and "L" for the low frequency setting. |
low_bound |
positive real number |
up_bound |
positive real number |
In this algorithm the preliminary estimate of alpha is found via using uniroot function. The latter is
given the lower and the upper bounds for alpha via low_bound and up_bound parameters. It is not possible to
pass 0 as the lower bound because there are numerical limitations on the alpha estimate, caused by the
length of the sample path and by numerical errors. p and p_prime must belong to the interval (0,1/2) (in the notation kept in rlfsm package)
The two powers cannot be equal.
Mazur S, Otryakhin D, Podolskij M (2020). “Estimation of the linear fractional stable motion.” Bernoulli, 26(1), 226–252. https://doi.org/10.3150/19-BEJ1124.
m<-45; M<-60; N<-2^10-M sigma<-0.3 p<-0.2; p_prime<-0.4 #### Continuous case lfsm<-path(N=N,m=m,M=M,alpha=1.8,H=0.8, sigma=sigma,freq='L',disable_X=FALSE,seed=3)$lfsm GenHighEstim(p=p,p_prime=p_prime,path=lfsm,freq="H") #### H-1/alpha<0 case lfsm<-path(N=N,m=m,M=M,alpha=0.8,H=0.8, sigma=sigma,freq='H',disable_X=FALSE,seed=3)$lfsm GenHighEstim(p=p,p_prime=p_prime,path=lfsm,freq="H")m<-45; M<-60; N<-2^10-M sigma<-0.3 p<-0.2; p_prime<-0.4 #### Continuous case lfsm<-path(N=N,m=m,M=M,alpha=1.8,H=0.8, sigma=sigma,freq='L',disable_X=FALSE,seed=3)$lfsm GenHighEstim(p=p,p_prime=p_prime,path=lfsm,freq="H") #### H-1/alpha<0 case lfsm<-path(N=N,m=m,M=M,alpha=0.8,H=0.8, sigma=sigma,freq='H',disable_X=FALSE,seed=3)$lfsm GenHighEstim(p=p,p_prime=p_prime,path=lfsm,freq="H")
General estimation procedure for low frequency case when 1/alpha is not a natural number.
GenLowEstim(t1, t2, p, path, freq = "L")GenLowEstim(t1, t2, p, path, freq = "L")
t1, t2
|
real number such that t2 > t1 > 0 |
p |
power |
path |
sample path of lfsm on which the inference is to be performed |
freq |
Frequency of the motion. It can take two values: "H" for high frequency and "L" for the low frequency setting. |
Mazur S, Otryakhin D, Podolskij M (2020). “Estimation of the linear fractional stable motion.” Bernoulli, 26(1), 226–252. https://doi.org/10.3150/19-BEJ1124.
m<-45; M<-60; N<-2^10-M sigma<-0.3 p<-0.3; k=3; t1=1; t2=2 #### Continuous case lfsm<-path(N=N,m=m,M=M,alpha=1.8,H=0.8, sigma=sigma,freq='L',disable_X=FALSE,seed=3)$lfsm GenLowEstim(t1=t1,t2=t2,p=p,path=lfsm,freq="L") #### H-1/alpha<0 case lfsm<-path(N=N,m=m,M=M,alpha=0.8,H=0.8, sigma=sigma,freq='L',disable_X=FALSE,seed=3)$lfsm GenLowEstim(t1=t1,t2=t2,p=p,path=lfsm,freq="L") #### The procedure works also for high frequency case lfsm<-path(N=N,m=m,M=M,alpha=1.8,H=0.8, sigma=sigma,freq='H',disable_X=FALSE,seed=3)$lfsm GenLowEstim(t1=t1,t2=t2,p=p,path=lfsm,freq="H")m<-45; M<-60; N<-2^10-M sigma<-0.3 p<-0.3; k=3; t1=1; t2=2 #### Continuous case lfsm<-path(N=N,m=m,M=M,alpha=1.8,H=0.8, sigma=sigma,freq='L',disable_X=FALSE,seed=3)$lfsm GenLowEstim(t1=t1,t2=t2,p=p,path=lfsm,freq="L") #### H-1/alpha<0 case lfsm<-path(N=N,m=m,M=M,alpha=0.8,H=0.8, sigma=sigma,freq='L',disable_X=FALSE,seed=3)$lfsm GenLowEstim(t1=t1,t2=t2,p=p,path=lfsm,freq="L") #### The procedure works also for high frequency case lfsm<-path(N=N,m=m,M=M,alpha=1.8,H=0.8, sigma=sigma,freq='H',disable_X=FALSE,seed=3)$lfsm GenLowEstim(t1=t1,t2=t2,p=p,path=lfsm,freq="H")
The statistic is defined as
H_hat(p, k, path)H_hat(p, k, path)
p |
power |
k |
increment order |
path |
sample path of lfsm on which the inference is to be performed |
Mazur S, Otryakhin D, Podolskij M (2020). “Estimation of the linear fractional stable motion.” Bernoulli, 26(1), 226–252. https://doi.org/10.3150/19-BEJ1124.
Function is given by
h_kr(k, r, x, H, alpha, l = 0)h_kr(k, r, x, H, alpha, l = 0)
k |
order of the increment, a natural number |
r |
difference step, a natural number |
x |
real number |
H |
Hurst parameter |
alpha |
self-similarity parameter of alpha stable random motion. |
l |
a value by which we shift x. Is used for computing function f_.+l and is passed to integrate function. |
Mazur S, Otryakhin D, Podolskij M (2020). “Estimation of the linear fractional stable motion.” Bernoulli, 26(1), 226–252. https://doi.org/10.3150/19-BEJ1124.
#### Plot h_kr #### s<-seq(0,10, by=0.01) h_val<-sapply(s,h_kr, k=5, r=1, H=0.3, alpha=1) plot(s,h_val)#### Plot h_kr #### s<-seq(0,10, by=0.01) h_val<-sapply(s,h_kr, k=5, r=1, H=0.3, alpha=1) plot(s,h_val)
Difference of the kth order. Defined as following:
Index i here is a coordinate in terms of point_num. Although R uses vector indexes that start from 1, increment has i varying from 0 to N, so that a vector has a length N+1. It is done in order to comply with the notation of the paper. This function doesn't allow for choosing frequency n. The frequency is determined by the path supplied, thus n equals to either the length of the path in high frequency setting or 1 in low frequency setting. increment() gives increments at certain point passed as i, which is a vector here. increments() computes high order increments for the whole sample path. The first function evaluates the formula above, while the second one uses structure diff(diff(...)) because the formula is slower at higher k.
increment(r, i, k, path) increments(k, r, path)increment(r, i, k, path) increments(k, r, path)
r |
difference step, a natural number |
i |
index of the point at which the increment is to be computed, a natural number. |
k |
order of the increment, a natural number |
path |
sample path for which a kth order increment is computed |
Mazur S, Otryakhin D, Podolskij M (2020). “Estimation of the linear fractional stable motion.” Bernoulli, 26(1), 226–252. https://doi.org/10.3150/19-BEJ1124.
m<-45; M<-60; N<-2^10-M alpha<-0.8; H<-0.8; sigma<-0.3 lfsm<-path(N=N,m=m,M=M,alpha=alpha,H=H, sigma=sigma,freq='L',disable_X=FALSE,seed=3)$lfsm tryCatch( increment(r=1,i=length(lfsm),k=length(lfsm)+100,path=lfsm), error=function(c) 'An error occures when k is larger then the length of the sample path') increment(r=3,i=50,k=3,path=lfsm) path=c(1,4,3,6,8,5,3,5,8,5,1,8,6) r=2; k=3 n <- length(path) - 1 DeltaX = increment(seq(r*k, n), path = path, k = k, r = r) DeltaX == increments(k=k,r=r,path)m<-45; M<-60; N<-2^10-M alpha<-0.8; H<-0.8; sigma<-0.3 lfsm<-path(N=N,m=m,M=M,alpha=alpha,H=H, sigma=sigma,freq='L',disable_X=FALSE,seed=3)$lfsm tryCatch( increment(r=1,i=length(lfsm),k=length(lfsm)+100,path=lfsm), error=function(c) 'An error occures when k is larger then the length of the sample path') increment(r=3,i=50,k=3,path=lfsm) path=c(1,4,3,6,8,5,3,5,8,5,1,8,6) r=2; k=3 n <- length(path) - 1 DeltaX = increment(seq(r*k, n), path = path, k = k, r = r) DeltaX == increments(k=k,r=r,path)
defined as for positive powers. When p is negative (-p is positive) the equality does not hold.
m_pk(k, p, alpha, H, sigma)m_pk(k, p, alpha, H, sigma)
k |
increment order |
p |
a positive number |
alpha |
self-similarity parameter of alpha stable random motion. |
H |
Hurst parameter |
sigma |
Scale parameter of lfsm |
The following identity is used for computations:
Mazur S, Otryakhin D, Podolskij M (2020). “Estimation of the linear fractional stable motion.” Bernoulli, 26(1), 226–252. https://doi.org/10.3150/19-BEJ1124.
The function is useful, for instance, when one needs to compute standard deviation of
estimator given a fixed set of parameters.
MCestimLFSM(Nmc, s, m, M, alpha, H, sigma, fr, Inference, ...)MCestimLFSM(Nmc, s, m, M, alpha, H, sigma, fr, Inference, ...)
Nmc |
Number of Monte Carlo repetitions |
s |
sequence of path lengths |
m |
discretization. A number of points between two nearby motion points |
M |
truncation parameter. A number of points at which the integral representing the definition of lfsm is calculated. So, after M points back we consider the rest of the integral to be 0. |
alpha |
self-similarity parameter of alpha stable random motion. |
H |
Hurst parameter |
sigma |
Scale parameter of lfsm |
fr |
frequency. Either "H" or "L" |
Inference |
statistical function to apply to sample paths |
... |
parameters to pass to Inference |
MCestimLFSM performs Monte-Carlo experiments to compute parameters according to procedure Inference. More specifically, for each element of s it generates Nmc lfsm sample paths with length equal to s[i], performs the statistical inference on each, obtaining the estimates, and then returns their different statistics. It is vital that the estimator returns a list of named parameters (one or several of 'sigma', 'alpha' and 'H'). MCestimLFSM uses the names to lookup the true parameter value and compute its bias.
For sample path generation MCestimLFSM uses a light-weight version of path, path_fast. In order to be applied, function Inference must accept argument 'path' as a sample path.
It returns a list containing the following components:
data |
a data frame, values of the estimates depending on path length s |
data_nor |
a data frame, normalized values of the estimates depending on path length s |
means, biases, sds
|
data frames: means, biases and standard deviations of the estimators depending on s |
Inference |
a function used to obtain estimates |
alpha, H, sigma
|
the parameters for which MCestimLFSM performs path generation |
freq |
frequency, either 'L' for low- or 'H' for high frequency |
#### Set of global parameters #### m<-25; M<-60 p<-.4; p_prime<-.2; k<-2 t1<-1; t2<-2 NmonteC<-5e1 S<-c(1e2,3e2) alpha<-1.8; H<-0.8; sigma<-0.3 # How to plot empirical density theor_3_1_H_clt<-MCestimLFSM(s=S,fr='H',Nmc=NmonteC, m=m,M=M,alpha=alpha,H=H, sigma=sigma,ContinEstim, t1=t1,t2=t2,p=p,k=k) l_plot<-Plot_dens(par_vec=c('sigma','alpha','H'), MC_data=theor_3_1_H_clt, Nnorm=1e7) # For MCestimLFSM() it is vital that the estimator returns a list of named parameters H_hat_f <- function(p,k,path) {hh<-H_hat(p,k,path); list(H=hh)} theor_3_1_H_clt<-MCestimLFSM(s=S,fr='H',Nmc=NmonteC, m=m,M=M,alpha=alpha,H=H, sigma=sigma,H_hat_f, p=p,k=k) # The estimator can return one, two or three of the parameters. est_1 <- function(path) list(H=1) theor_3_1_H_clt<-MCestimLFSM(s=S,fr='H',Nmc=NmonteC, m=m,M=M,alpha=alpha,H=H, sigma=sigma,est_1) est_2 <- function(path) list(H=0.8, alpha=1.5) theor_3_1_H_clt<-MCestimLFSM(s=S,fr='H',Nmc=NmonteC, m=m,M=M,alpha=alpha,H=H, sigma=sigma,est_2) est_3 <- function(path) list(sigma=5, H=0.8, alpha=1.5) theor_3_1_H_clt<-MCestimLFSM(s=S,fr='H',Nmc=NmonteC, m=m,M=M,alpha=alpha,H=H, sigma=sigma,est_3)#### Set of global parameters #### m<-25; M<-60 p<-.4; p_prime<-.2; k<-2 t1<-1; t2<-2 NmonteC<-5e1 S<-c(1e2,3e2) alpha<-1.8; H<-0.8; sigma<-0.3 # How to plot empirical density theor_3_1_H_clt<-MCestimLFSM(s=S,fr='H',Nmc=NmonteC, m=m,M=M,alpha=alpha,H=H, sigma=sigma,ContinEstim, t1=t1,t2=t2,p=p,k=k) l_plot<-Plot_dens(par_vec=c('sigma','alpha','H'), MC_data=theor_3_1_H_clt, Nnorm=1e7) # For MCestimLFSM() it is vital that the estimator returns a list of named parameters H_hat_f <- function(p,k,path) {hh<-H_hat(p,k,path); list(H=hh)} theor_3_1_H_clt<-MCestimLFSM(s=S,fr='H',Nmc=NmonteC, m=m,M=M,alpha=alpha,H=H, sigma=sigma,H_hat_f, p=p,k=k) # The estimator can return one, two or three of the parameters. est_1 <- function(path) list(H=1) theor_3_1_H_clt<-MCestimLFSM(s=S,fr='H',Nmc=NmonteC, m=m,M=M,alpha=alpha,H=H, sigma=sigma,est_1) est_2 <- function(path) list(H=0.8, alpha=1.5) theor_3_1_H_clt<-MCestimLFSM(s=S,fr='H',Nmc=NmonteC, m=m,M=M,alpha=alpha,H=H, sigma=sigma,est_2) est_3 <- function(path) list(sigma=5, H=0.8, alpha=1.5) theor_3_1_H_clt<-MCestimLFSM(s=S,fr='H',Nmc=NmonteC, m=m,M=M,alpha=alpha,H=H, sigma=sigma,est_3)
Alpha-norm of an arbitrary function
Norm_alpha(fun, alpha, ...)Norm_alpha(fun, alpha, ...)
fun |
a function to compute a norm |
alpha |
self-similarity parameter of alpha stable random motion. |
... |
a set of parameters to pass to integrate |
fun must accept a vector of values for evaluation. See ?integrate for further details. Most problems with this function appear because of rather high precision. Try to tune rel.tol parameter first.
Mazur S, Otryakhin D, Podolskij M (2020). “Estimation of the linear fractional stable motion.” Bernoulli, 26(1), 226–252. https://doi.org/10.3150/19-BEJ1124.
Norm_alpha(h_kr,alpha=1.8,k=2,r=1,H=0.8,l=4)Norm_alpha(h_kr,alpha=1.8,k=2,r=1,H=0.8,l=4)
The function creates a 1-dimensional LFSM sample path using the numerical algorithm from the paper by Otryakhin and Mazur. The theoretical foundation of the method comes from the article by Stoev and Taqqu. Linear fractional stable motion is defined as
path( N = NULL, m, M, alpha, H, sigma, freq, disable_X = FALSE, levy_increments = NULL, seed = NULL )path( N = NULL, m, M, alpha, H, sigma, freq, disable_X = FALSE, levy_increments = NULL, seed = NULL )
N |
a number of points of the lfsm. |
m |
discretization. A number of points between two nearby motion points |
M |
truncation parameter. A number of points at which the integral representing the definition of lfsm is calculated. So, after M points back we consider the rest of the integral to be 0. |
alpha |
self-similarity parameter of alpha stable random motion. |
H |
Hurst parameter |
sigma |
Scale parameter of lfsm |
freq |
Frequency of the motion. It can take two values: "H" for high frequency and "L" for the low frequency setting. |
disable_X |
is needed to disable computation of X. The default value is FALSE. When it is TRUE, only a levy motion is returned, which in turn reduces the computation time. The feature is particularly useful for reproducibility when combined with seeding. |
levy_increments |
increments of Levy motion underlying the lfsm. |
seed |
this parameter performs seeding of path generator |
It returns a list containing the motion, the underlying Levy motion, the point number of the motions from 0 to N and the corresponding coordinate (which depends on the frequency), the parameters that were used to generate the lfsm, and the predefined frequency.
Mazur S, Otryakhin D (2020). “Linear Fractional Stable Motion with the rlfsm R Package.” The R Journal, 12(1), 386–405. doi:10.32614/RJ-2020-008.
Stoev S, Taqqu MS (2004). “Simulation methods for linear fractional stable motion and FARIMA using the Fast Fourier Transform.” Fractals, 95(1), 95-121. https://doi.org/10.1142/S0218348X04002379.
paths simulates a number of lfsm sample paths.
# Path generation m<-256; M<-600; N<-2^10-M alpha<-1.8; H<-0.8; sigma<-0.3 seed=2 List<-path(N=N,m=m,M=M,alpha=alpha,H=H, sigma=sigma,freq='L',disable_X=FALSE,seed=3) # Normalized paths Norm_lfsm<-List[['lfsm']]/max(abs(List[['lfsm']])) Norm_oLm<-List[['levy_motion']]/max(abs(List[['levy_motion']])) # Visualization of the paths plot(Norm_lfsm, col=2, type="l", ylab="coordinate") lines(Norm_oLm, col=3) leg.txt <- c("lfsm", "oLm") legend("topright",legend = leg.txt, col =c(2,3), pch=1) # Creating Levy motion levyIncrems<-path(N=N, m=m, M=M, alpha, H, sigma, freq='L', disable_X=TRUE, levy_increments=NULL, seed=seed) # Creating lfsm based on the levy motion lfsm_full<-path(m=m, M=M, alpha=alpha, H=H, sigma=sigma, freq='L', disable_X=FALSE, levy_increments=levyIncrems$levy_increments, seed=seed) sum(levyIncrems$levy_increments== lfsm_full$levy_increments)==length(lfsm_full$levy_increments)# Path generation m<-256; M<-600; N<-2^10-M alpha<-1.8; H<-0.8; sigma<-0.3 seed=2 List<-path(N=N,m=m,M=M,alpha=alpha,H=H, sigma=sigma,freq='L',disable_X=FALSE,seed=3) # Normalized paths Norm_lfsm<-List[['lfsm']]/max(abs(List[['lfsm']])) Norm_oLm<-List[['levy_motion']]/max(abs(List[['levy_motion']])) # Visualization of the paths plot(Norm_lfsm, col=2, type="l", ylab="coordinate") lines(Norm_oLm, col=3) leg.txt <- c("lfsm", "oLm") legend("topright",legend = leg.txt, col =c(2,3), pch=1) # Creating Levy motion levyIncrems<-path(N=N, m=m, M=M, alpha, H, sigma, freq='L', disable_X=TRUE, levy_increments=NULL, seed=seed) # Creating lfsm based on the levy motion lfsm_full<-path(m=m, M=M, alpha=alpha, H=H, sigma=sigma, freq='L', disable_X=FALSE, levy_increments=levyIncrems$levy_increments, seed=seed) sum(levyIncrems$levy_increments== lfsm_full$levy_increments)==length(lfsm_full$levy_increments)
The function takes a list of parameters (alpha, H) and uses expand.grid to obtain all possible combinations of them.
Based on each combination, the function simulates an lfsm sample path. It is meant to be used in conjunction with function Plot_list_paths.
Path_array(N, m, M, l, sigma)Path_array(N, m, M, l, sigma)
N |
a number of points of the lfsm. |
m |
discretization. A number of points between two nearby motion points |
M |
truncation parameter. A number of points at which the integral representing the definition of lfsm is calculated. So, after M points back we consider the rest of the integral to be 0. |
l |
a list of parameters to expand |
sigma |
Scale parameter of lfsm |
The returned value is a data frame containing paths and the corresponding values of alpha, H and frequency.
l=list(H=c(0.2,0.8),alpha=c(1,1.8), freq="H") arr<-Path_array(N=300,m=30,M=100,l=l,sigma=0.3) str(arr) head(arr)l=list(H=c(0.2,0.8),alpha=c(1,1.8), freq="H") arr<-Path_array(N=300,m=30,M=100,l=l,sigma=0.3) str(arr) head(arr)
It is essentially a wrapper for path generator, which exploits the latest to create a matrix with paths in its columns.
paths(N_var, parallel, seed_list = rep(x = NULL, times = N_var), ...)paths(N_var, parallel, seed_list = rep(x = NULL, times = N_var), ...)
N_var |
number of lfsm paths to generate |
parallel |
a TRUE/FALSE flag which determines if the paths will be created in parallel or sequentially |
seed_list |
a numerical vector of seeds to pass to |
... |
arguments to pass to path |
m<-45; M<-60; N<-2^10-M alpha<-1.8; H<-0.8; sigma<-0.3 freq='L' r=1; k=2; p=0.4 Y<-paths(N_var=10,parallel=TRUE,N=N,m=m,M=M, alpha=alpha,H=H,sigma=sigma,freq='L', disable_X=FALSE,levy_increments=NULL) Hs<-apply(Y,MARGIN=2,H_hat,p=p,k=k) hist(Hs)m<-45; M<-60; N<-2^10-M alpha<-1.8; H<-0.8; sigma<-0.3 freq='L' r=1; k=2; p=0.4 Y<-paths(N_var=10,parallel=TRUE,N=N,m=m,M=M, alpha=alpha,H=H,sigma=sigma,freq='L', disable_X=FALSE,levy_increments=NULL) Hs<-apply(Y,MARGIN=2,H_hat,p=p,k=k) hist(Hs)
Defined as
,
where
phi(t, k, path, H, freq)phi(t, k, path, H, freq)
t |
positive real number |
k |
increment order |
path |
sample path of lfsm on which the inference is to be performed |
H |
Hurst parameter |
freq |
Frequency of the motion. It can take two values: "H" for high frequency and "L" for the low frequency setting. |
Hurst parameter is required only in high frequency case. In the low frequency, there is no need to assign H a value because it will not be evaluated.
Mazur S, Otryakhin D, Podolskij M (2020). “Estimation of the linear fractional stable motion.” Bernoulli, 26(1), 226–252. https://doi.org/10.3150/19-BEJ1124.
A function from a general estimation procedure which is defined as m^p_-p'_k /m^p'_-p_k, originally proposed in [13].
phi_of_alpha(p, p_prime, alpha)phi_of_alpha(p, p_prime, alpha)
p |
power |
p_prime |
power |
alpha |
self-similarity parameter of alpha stable random motion. |
Mazur S, Otryakhin D, Podolskij M (2020). “Estimation of the linear fractional stable motion.” Bernoulli, 26(1), 226–252. https://doi.org/10.3150/19-BEJ1124.
Plots the densities of the parameters (alpha,H,sigma) estimated in Monte-Carlo experiment.
Works in conjunction with MCestimLFSM function.
Plot_dens(par_vec = c("alpha", "H", "sigma"), MC_data, Nnorm = 1e+07)Plot_dens(par_vec = c("alpha", "H", "sigma"), MC_data, Nnorm = 1e+07)
par_vec |
vector of parameters which are to be plotted |
MC_data |
a list created by |
Nnorm |
number of point sampled from standard normal distribution |
Plot_vb to plot variance- and bias dependencies on n.
m<-45; M<-60 p<-.4; p_prime<-.2 t1<-1; t2<-2; k<-2 NmonteC<-5e2 S<-c(1e3,1e4) alpha<-.8; H<-0.8; sigma<-0.3 theor_4_1_clt_new<-MCestimLFSM(s=S,fr='L',Nmc=NmonteC, m=m,M=M, alpha=alpha,H=H,sigma=sigma, GenLowEstim,t1=t1,t2=t2,p=p) l_plot<-Plot_dens(par_vec=c('sigma','alpha','H'), MC_data=theor_4_1_clt_new, Nnorm=1e7) l_plotm<-45; M<-60 p<-.4; p_prime<-.2 t1<-1; t2<-2; k<-2 NmonteC<-5e2 S<-c(1e3,1e4) alpha<-.8; H<-0.8; sigma<-0.3 theor_4_1_clt_new<-MCestimLFSM(s=S,fr='L',Nmc=NmonteC, m=m,M=M, alpha=alpha,H=H,sigma=sigma, GenLowEstim,t1=t1,t2=t2,p=p) l_plot<-Plot_dens(par_vec=c('sigma','alpha','H'), MC_data=theor_4_1_clt_new, Nnorm=1e7) l_plot
Rendering of path lattice
Plot_list_paths(arr)Plot_list_paths(arr)
arr |
a data frame produced by |
l=list(H=c(0.2,0.8),alpha=c(1,1.8), freq="H") arr<-Path_array(N=300,m=30,M=100,l=l,sigma=0.3) p<-Plot_list_paths(arr) pl=list(H=c(0.2,0.8),alpha=c(1,1.8), freq="H") arr<-Path_array(N=300,m=30,M=100,l=l,sigma=0.3) p<-Plot_list_paths(arr) p
MCestimLFSM function.A function to plot variance- and bias dependencies of estimators on the lengths of sample paths. Works in conjunction with MCestimLFSM function.
Plot_vb(data)Plot_vb(data)
data |
a list created by |
The function returns a ggplot2 graph.
# Light weight computaions m<-25; M<-50 alpha<-1.8; H<-0.8; sigma<-0.3 S<-c(1:3)*1e2 p<-.4; p_prime<-.2; t1<-1; t2<-2 k<-2; NmonteC<-50 # Here is the continuous H-1/alpha inference procedure theor_3_1_H_clt<-MCestimLFSM(s=S,fr='H',Nmc=NmonteC, m=m,M=M,alpha=alpha,H=H, sigma=sigma,ContinEstim, t1=t1,t2=t2,p=p,k=k) Plot_vb(theor_3_1_H_clt) # More demanding example (it is better to use multicore setup) # General low frequency inference m<-45; M<-60 alpha<-0.8; H<-0.8; sigma<-0.3 S<-c(1:15)*1e2 p<-.4; t1<-1; t2<-2 NmonteC<-50 # Here is the continuous H-1/alpha inference procedure theor_4_1_H_clt<-MCestimLFSM(s=S,fr='H',Nmc=NmonteC, m=m,M=M,alpha=alpha,H=H, sigma=sigma,GenLowEstim, t1=t1,t2=t2,p=p) Plot_vb(theor_4_1_H_clt)# Light weight computaions m<-25; M<-50 alpha<-1.8; H<-0.8; sigma<-0.3 S<-c(1:3)*1e2 p<-.4; p_prime<-.2; t1<-1; t2<-2 k<-2; NmonteC<-50 # Here is the continuous H-1/alpha inference procedure theor_3_1_H_clt<-MCestimLFSM(s=S,fr='H',Nmc=NmonteC, m=m,M=M,alpha=alpha,H=H, sigma=sigma,ContinEstim, t1=t1,t2=t2,p=p,k=k) Plot_vb(theor_3_1_H_clt) # More demanding example (it is better to use multicore setup) # General low frequency inference m<-45; M<-60 alpha<-0.8; H<-0.8; sigma<-0.3 S<-c(1:15)*1e2 p<-.4; t1<-1; t2<-2 NmonteC<-50 # Here is the continuous H-1/alpha inference procedure theor_4_1_H_clt<-MCestimLFSM(s=S,fr='H',Nmc=NmonteC, m=m,M=M,alpha=alpha,H=H, sigma=sigma,GenLowEstim, t1=t1,t2=t2,p=p) Plot_vb(theor_4_1_H_clt)
Defined as
R_hl(p, k, path)R_hl(p, k, path)
p |
power |
k |
increment order |
path |
sample path of lfsm on which the inference is to be performed |
The computation procedure for high- and low frequency cases is the same, since there is no way to control frequency given a sample path.
Mazur S, Otryakhin D, Podolskij M (2020). “Estimation of the linear fractional stable motion.” Bernoulli, 26(1), 226–252. https://doi.org/10.3150/19-BEJ1124.
m<-45; M<-60; N<-2^10-M alpha<-0.8; H<-0.8; sigma<-0.3 p<-0.3; k=3 lfsm<-path(N=N,m=m,M=M,alpha=alpha,H=H, sigma=sigma,freq='L',disable_X=FALSE,seed=3)$lfsm R_hl(p=p,k=k,path=lfsm)m<-45; M<-60; N<-2^10-M alpha<-0.8; H<-0.8; sigma<-0.3 p<-0.3; k=3 lfsm<-path(N=N,m=m,M=M,alpha=alpha,H=H, sigma=sigma,freq='L',disable_X=FALSE,seed=3)$lfsm R_hl(p=p,k=k,path=lfsm)
Retrieve statistics(bias, variance) of estimators based on a set of paths
Retrieve_stats(paths, true_val, Est, ...)Retrieve_stats(paths, true_val, Est, ...)
paths |
real-valued matrix representing sample paths of the stochastic process being studied |
true_val |
true value of the estimated parameter |
Est |
estimator (i.e. H_hat) |
... |
parameters to pass to Est |
m<-45; M<-60; N<-2^10-M alpha<-1.8; H<-0.8; sigma<-0.3 freq='L';t1=1; t2=2 r=1; k=2; p=0.4 Y<-paths(N_var=10,parallel=TRUE,N=N,m=m,M=M, alpha=alpha,H=H,sigma=sigma,freq='L', disable_X=FALSE,levy_increments=NULL) Retrieve_stats(paths=Y,true_val=sigma,Est=sigma_hat,t1=t1,k=2,alpha=alpha,H=H,freq="L")m<-45; M<-60; N<-2^10-M alpha<-1.8; H<-0.8; sigma<-0.3 freq='L';t1=1; t2=2 r=1; k=2; p=0.4 Y<-paths(N_var=10,parallel=TRUE,N=N,m=m,M=M, alpha=alpha,H=H,sigma=sigma,freq='L', disable_X=FALSE,levy_increments=NULL) Retrieve_stats(paths=Y,true_val=sigma,Est=sigma_hat,t1=t1,k=2,alpha=alpha,H=H,freq="L")
Statistic of the form
sf(path, f, k, r, H, freq, ...)sf(path, f, k, r, H, freq, ...)
path |
sample path for which the statistic is to be calculated. |
f |
function applied to high order increments. |
k |
order of the increments. |
r |
step of high order increments. |
H |
Hurst parameter. |
freq |
frequency. |
... |
parameters to pass to function f |
Hurst parameter is required only in high frequency case. In the low frequency, there is no need to assign H a value because it will not be evaluated.
Mazur S, Otryakhin D, Podolskij M (2020). “Estimation of the linear fractional stable motion.” Bernoulli, 26(1), 226–252. https://doi.org/10.3150/19-BEJ1124.
phi computes V statistic with f(.)=cos(t.)
m<-45; M<-60; N<-2^10-M alpha<-1.8; H<-0.8; sigma<-0.3 freq='L' r=1; k=2; p=0.4 S<-(1:20)*100 path_lfsm<-function(...){ List<-path(...) List$lfsm } Pths<-lapply(X=S,FUN=path_lfsm, m=m, M=M, alpha=alpha, sigma=sigma, H=H, freq=freq, disable_X = FALSE, levy_increments = NULL, seed = NULL) f_phi<-function(t,x) cos(t*x) f_pow<-function(x,p) (abs(x))^p V_cos<-sapply(Pths,FUN=sf,f=f_phi,k=k,r=r,H=H,freq=freq,t=1) ex<-exp(-(abs(sigma*Norm_alpha(h_kr,alpha=alpha,k=k,r=r,H=H,l=0)$result)^alpha)) # Illustration of the law of large numbers for phi: plot(y=V_cos, x=S, ylim = c(0,max(V_cos)+0.1)) abline(h=ex, col='brown') # Illustration of the law of large numbers for power functions: Mpk<-m_pk(k=k, p=p, alpha=alpha, H=H, sigma=sigma) sf_mod<-function(Xpath,...) { Path<-unlist(Xpath) sf(path=Path,...) } V_pow<-sapply(Pths,FUN=sf_mod,f=f_pow,k=k,r=r,H=H,freq=freq,p=p) plot(y=V_pow, x=S, ylim = c(0,max(V_pow)+0.1)) abline(h=Mpk, col='brown')m<-45; M<-60; N<-2^10-M alpha<-1.8; H<-0.8; sigma<-0.3 freq='L' r=1; k=2; p=0.4 S<-(1:20)*100 path_lfsm<-function(...){ List<-path(...) List$lfsm } Pths<-lapply(X=S,FUN=path_lfsm, m=m, M=M, alpha=alpha, sigma=sigma, H=H, freq=freq, disable_X = FALSE, levy_increments = NULL, seed = NULL) f_phi<-function(t,x) cos(t*x) f_pow<-function(x,p) (abs(x))^p V_cos<-sapply(Pths,FUN=sf,f=f_phi,k=k,r=r,H=H,freq=freq,t=1) ex<-exp(-(abs(sigma*Norm_alpha(h_kr,alpha=alpha,k=k,r=r,H=H,l=0)$result)^alpha)) # Illustration of the law of large numbers for phi: plot(y=V_cos, x=S, ylim = c(0,max(V_cos)+0.1)) abline(h=ex, col='brown') # Illustration of the law of large numbers for power functions: Mpk<-m_pk(k=k, p=p, alpha=alpha, H=H, sigma=sigma) sf_mod<-function(Xpath,...) { Path<-unlist(Xpath) sf(path=Path,...) } V_pow<-sapply(Pths,FUN=sf_mod,f=f_pow,k=k,r=r,H=H,freq=freq,p=p) plot(y=V_pow, x=S, ylim = c(0,max(V_pow)+0.1)) abline(h=Mpk, col='brown')
Statistical estimator for sigma
sigma_hat(t1, k, path, alpha, H, freq)sigma_hat(t1, k, path, alpha, H, freq)
t1 |
real number such that t1 > 0 |
k |
increment order |
path |
sample path of lfsm on which the inference is to be performed |
alpha |
self-similarity parameter of alpha stable random motion. |
H |
Hurst parameter |
freq |
Frequency of the motion. It can take two values: "H" for high frequency and "L" for the low frequency setting. |
m<-45; M<-60; N<-2^14-M alpha<-1.8; H<-0.8; sigma<-0.3 freq='H' r=1; k=2; p=0.4; t1=1; t2=2 # Reproducing the work of ContinEstim # in high frequency case lfsm<-path(N=N,m=m,M=M,alpha=alpha,H=H, sigma=sigma,freq='L',disable_X=FALSE,seed=1)$lfsm H_est<-H_hat(p=p,k=k,path=lfsm) H_est alpha_est<-alpha_hat(t1=t1,t2=t2,k=k,path=lfsm,H=H_est,freq=freq) alpha_est sigma_est<-tryCatch( sigma_hat(t1=t1,k=k,path=lfsm, alpha=alpha_est,H=H_est,freq=freq), error=function(c) 'Impossible to compute sigma_est') sigma_estm<-45; M<-60; N<-2^14-M alpha<-1.8; H<-0.8; sigma<-0.3 freq='H' r=1; k=2; p=0.4; t1=1; t2=2 # Reproducing the work of ContinEstim # in high frequency case lfsm<-path(N=N,m=m,M=M,alpha=alpha,H=H, sigma=sigma,freq='L',disable_X=FALSE,seed=1)$lfsm H_est<-H_hat(p=p,k=k,path=lfsm) H_est alpha_est<-alpha_hat(t1=t1,t2=t2,k=k,path=lfsm,H=H_est,freq=freq) alpha_est sigma_est<-tryCatch( sigma_hat(t1=t1,k=k,path=lfsm, alpha=alpha_est,H=H_est,freq=freq), error=function(c) 'Impossible to compute sigma_est') sigma_est
Function of the form
theta(p, alpha, sigma, g, h)theta(p, alpha, sigma, g, h)
p |
power, real number from (-1,1) |
alpha |
self-similarity parameter of alpha stable random motion. |
sigma |
Scale parameter of lfsm |
g, h
|
functions |
Mazur S, Otryakhin D, Podolskij M (2020). “Estimation of the linear fractional stable motion.” Bernoulli, 26(1), 226–252. https://doi.org/10.3150/19-BEJ1124.
alpha norm of u*g
U_g(g, u, ...)U_g(g, u, ...)
g |
function |
u |
real number |
... |
additional parameters to pass to Norm_alpha |
g<-function(x) exp(-x^2) g<-function(x) exp(-abs(x)) U_g(g=g,u=4,alpha=1.7)g<-function(x) exp(-x^2) g<-function(x) exp(-abs(x)) U_g(g=g,u=4,alpha=1.7)
alpha-norm of u*g + v*h.
U_gh(g, h, u, v, ...)U_gh(g, h, u, v, ...)
g, h
|
functions |
v, u
|
real numbers |
... |
additional parameters to pass to Norm_alpha |
g<-function(x) exp(-x^2) h<-function(x) exp(-abs(x)) U_gh(g=g, h=h, u=4, v=3, alpha=1.7)g<-function(x) exp(-x^2) h<-function(x) exp(-abs(x)) U_gh(g=g, h=h, u=4, v=3, alpha=1.7)
It is used when random variables do not have finite second moments, and thus, the covariance matrix is not defined.
For and with . Then the measure of dependence is given by via
U_ghuv(alpha, sigma, g, h, u, v, ...)U_ghuv(alpha, sigma, g, h, u, v, ...)
alpha |
self-similarity parameter of alpha stable random motion. |
sigma |
Scale parameter of lfsm |
g, h
|
functions |
v, u
|
real numbers |
... |
additional parameters to pass to U_gh and U_g |
g<-function(x) exp(-x^2) h<-function(x) exp(-abs(x)) U_ghuv(alpha=1.5, sigma=1, g=g, h=h, u=10, v=15, rel.tol = .Machine$double.eps^0.25, abs.tol=1e-11)g<-function(x) exp(-x^2) h<-function(x) exp(-abs(x)) U_ghuv(alpha=1.5, sigma=1, g=g, h=h, u=10, v=15, rel.tol = .Machine$double.eps^0.25, abs.tol=1e-11)