LongMemory.jl: Generating, Estimating, and Forecasting Long Memory Models in Julia

Illustrative example

Author

J. Eduardo Vera-Valdés, eduardo@math.aau.dk

Abstract

LongMemory.jl is a package for time series long memory modelling in Julia. The package provides functions for generating long memory, estimating the parameters of the models, and simulation. Generating methods include fractional differencing, stochastic error duration, and cross-sectional aggregation. Estimators include those inspired by the log-periodogram regression and parametric ones. This notebook provides an illustrative example on the use of the package.

Citation

If you use this package in your research, please cite it as:

Vera-Valdés, J.E. (2024). “LongMemory.jl: Generating, Estimating, and Forecasting Long Memory Models in Julia”. arXiv 2401.14077. https://arxiv.org/abs/2401.14077

@article{VERAVALDES2024a,
title = {LongMemory.jl: Generating, Estimating, and Forecasting Long Memory Models in Julia},
year = {2024},
author = {Vera-Valdés, J.E.},
journal = {arXiv preprint arXiv:2401.14077},
url = {https://arxiv.org/abs/2401.14077}
}

Nile River minima data

Plotting the data, its autocorrelation function and spectral density

Loading the packages.

using LongMemory, Plots, DataFrames, CSV

Data loading, and plot.

NileMin = NileData()
theme(:ggplot2)
p1 = plot( NileMin.Year , NileMin.NileMin, xlabel="Year" , ylabel = "Nile River minima" , legend = false )

Plotting the autocorrelation function using the autocorrelation_plot() function from LongMemory.jl.

p2 = autocorrelation_plot(NileMin.NileMin, 50)

Plotting the spectral density with the periodogram() function.

p3 = periodogram_plot(NileMin.NileMin, slope = true)

Combining the plots.

l = @layout [a; b c]
theme(:ggplot2)
plot(p1, p2, p3, layout = l, size = (700, 500) )

All the above plus the log-variance plot can be done in a single plot.

NileDataPlot()

Long memory generation

Fractional differencing generation and plotting

Loading the package for random number generation, setting the seed, and generating the fractional differenced series using the fi_gen() function.

using Random
Random.seed!(1234)

dx = fi_gen(500, 0.3)
500-element Vector{Float64}:
  0.9706563288552146
 -0.6880215128786354
  0.7973733442603006
  0.19192066674807962
 -0.46144511739243677
 -1.514799781765257
  2.2501256016290645
  2.0467649947720434
  1.5280686228970994
 -0.1269305253098677
  1.0031743811925413
  1.5910264204337021
  1.5990097697535517
  ⋮
  0.16749590571039433
 -1.1000791653801492
 -0.2066614957376673
 -0.9301341055606994
 -1.3332191794608566
 -0.4967820828987304
 -1.0341442518867014
 -1.0834385445422616
 -1.0557333271693792
 -0.08535062227423523
 -1.783593060065521
  0.11117399279145923

Plotting the fractional differenced series, and the sample autocorrelation function. The theoretical autocorrelation function is also plotted by using the fi_var_vals() function.

p1 = plot( dx , label = "Fractionally differenced data" )
p2 = plot( 0:50 , autocorrelation( dx , 51 ) , label = "Sample autocorrelation function", line = :stem , marker = :circle)
plot!( 0:50, fi_cor_vals( 51, 0.3 ), linewidth = 2, line = :dash, label = "Theoretical autocorrelation function")
l = @layout [a b]
theme(:ggplot2)
plot(p1, p2, layout = l, size = (700, 250) )

Cross-sectional aggregation generation and plotting

Generating the cross-sectional aggregation series using the csa_gen() function.

Multiple dispatch is used to generate the series. The finite approximation is used if we define the first two arguments as the sample size and the number of series. The asymptotic version is used if the number of series to use is not passed to the function. The final two arguments are the parameters of the Beta function.

csa_fin = csa_gen(1000,1000,1.3,1.5)
csa_asym = csa_gen(1000,1.3,1.5)
1000-element Vector{Float64}:
 -0.4367168712381375
 -0.27487348832431535
 -1.4597720231887006
 -1.9543596033122075
 -1.7243709928522584
 -2.4428086493018686
 -0.9191629868753005
  0.7415293874324638
  1.134795740715502
  0.6279072264178107
  0.41369790962276587
  0.08485016809582824
 -0.6521029064165207
  ⋮
  4.478681154035048
  3.081813222701477
  3.06039500311814
  3.1571288452388684
  1.9223537754308386
  1.9141232652309057
  2.469011881695831
  3.0677520738078097
  2.6416344731566843
  2.089792892141115
  1.6435055320286303
  1.2060575790942034

Plotting the autocorrelation function of the finite and asymptotic series. The theoretical autocorrelation function is also plotted by using the csa_var_vals() function.

plot( 0:50 , autocorrelation( csa_fin , 51 ) , label = "Sample autocorrelation function, finite approximation", line = :stem , marker = :circle)
plot!( 0:50 , autocorrelation( csa_asym , 51 ) , label = "Sample autocorrelation function, asymptotic model", line = :stem , marker = :circle, color = :black)
plot!( 0:50, csa_cor_vals( 51, 1.3, 1.5 ), linewidth = 3, line = :dash, label = "Theoretical autocorrelation function", color = :red)
theme(:ggplot2)
plot!(size = (500, 200) )

Stochastic duration shocks generation and plotting

Data generation and plotting. We use the sds_gen() function to generate the data, and the LMPlot() function to plot it along with three diagnostic plots.

LMPlot(sds_gen(1000,0.45), name = "Stochastic duration shock")
theme(:ggplot2)
plot!( size = (700, 500) )

Long memory estimation

Classic estimators

Variance plot in Nile River data

The variance plot is computed by the log_variance_plot() function. The number of subsamples to use is passed to the function as the m argument. The slope and slope2 arguments are used to plot the theoretical and estimated slopes, respectively.

log_variance_plot( NileData().NileMin; m=300, slope = true, slope2 = true )
plot!(size = (500, 300) )

Computing the long memory parameter by the log-variance plot.

log_variance_est( NileData().NileMin; m=300 )
0.15619757726603672

Rescaled range analysis in Nile River data

Estimation by the rescaled range, \(R/S\). The number of subsamples to use is passed to the function as the k argument.

rescaled_range_est( NileData().NileMin; k = 300 )
0.4254606013817649

Rescaled range plot. The number of subsamples to use is passed to the function as the k argument. The slope and slope2 arguments are used to plot the theoretical and estimated slopes, respectively.

rescaled_range_plot( NileData().NileMin; k = 300, slope = true, slope2 = true )

Semiparametric estimators

Log-periodogram regression in Nile River data

GPH and bias-reduced estimator with one polynomial term. The number of polynomial terms is passed to the function as the br argument. If not passed, the function uses the default value of 0, recovering the GPH estimator.

(gph_est( NileData().NileMin ),  gph_est( NileData().NileMin ; br = 1))
(0.37449410505423664, 0.39745526593583125)

Variance estimates using the data.

(gph_est_variance( NileData().NileMin ),  gph_est_variance( NileData().NileMin ; br = 1))
(0.002272008379624622, 0.0051120188541553995)

Variance estimates using multiple dispatch on the sample size.

(gph_est_variance( length( NileData().NileMin ) ),  gph_est_variance( length( NileData().NileMin ) ; br = 1))
(0.002272008379624622, 0.0051120188541553995)

(Exact) Local Whittle estimator in Nile River data

The local Whittle and exact local Whittle estimators are computed by the whittle_est() and exact_whittle_est() functions, respectively.

(whittle_est( NileData().NileMin ),  exact_whittle_est( NileData().NileMin ))
(0.37635955766433826, 0.4088495239569418)

Variance estimates using the data.

(whittle_est_variance( NileData().NileMin ),  exact_whittle_est_variance( NileData().NileMin ))
(0.0013812154696132596, 0.0013812154696132596)

Parametric estimators

Fractional differencing MLE in Nile River data

The fractional differencing MLE is computed by the fi_mle_est() function. The function returns the estimated fractional differencing parameter and the estimated standard error in a tuple.

dmle, sigfi = fi_mle_est( NileData().NileMin )
(0.3925714993964694, 69.95632676539786)

CSA MLE in Nile River data

The MLE for the cross-sectionally aggregated process is computed by the csa_mle_est() function. The function returns the estimated parameters of the Beta distribution, and the estimated standard error in a tuple.

pmle, qmle, sigcsa = csa_mle_est( NileData().NileMin )
(1.0000000537585416, 2.4475036737227627, 106.79292973709306)

Custom-HAR in Nile River data

A custom HAR model can be estimated by the function har_est(). The lags for the custom HAR model are passed as an array to the function as the m argument. The function returns the estimated parameters of the model, and the estimated standard error in a tuple.

beta_har, sighar = har_est( NileData().NileMin; m = [1, 7] )
([254.23541690816745, 0.40096895301134294, 0.377482428389992], 69.6124509836161)

Long memory forecasting

Forecasting the Nile River minima data

Fractional differencing forecasting

The fi_forecast_() function returns the forecasted values, while the fi_forecast_plot() functions returns the plot. First argument is the time series to forecast. The second and third arguments are used to define the number of steps ahead to forecast and the long memory parameter. The optional fourth argument is used to define the standard deviation to plot the confidence bands.

p1 = fi_forecast_plot( NileData().NileMin , 30, dmle, sigfi)

Cross-sectional aggregation forecasting

The csa_forecast_() function returns the forecasted values, while the csa_forecast_plot() functions returns the plot. The arguments are analogous to the fractional differencing plot function, noting that two Beta parameters are needed.

p2 = csa_forecast_plot( NileData().NileMin , 30, pmle, qmle, sigcsa)

Custom-HAR forecasting

The har_forecast_() function returns the forecasted values, while the har_forecast_plot() functions returns the plot. Nonetheless, the arguments are not analogous to the fractional differencing and cross-sectional aggregation plot functions, given that the custom-HAR model is estimated as part of the forecast.

Hence, the first argument is the time series to forecast and the second argument is used to define the number of steps ahead to forecast as before. The third optional argument, flag, is used to determine if the the confidence bands should be included in the plot. The fourth optional argument, m, is used to define the number of lags to use in the custom-HAR model. The optional fifth argument is used to define the standard deviation to plot the confidence bands.

p3 = har_forecast_plot( NileData().NileMin , 30; flag=true, m=[1, 7])

Comparing the forecasts

We generate a plot with the three forecasts using the @layout macro from the Plots.jl package.

l = @layout [a b c]
theme(:ggplot2)
plot(p1, p2, p3, layout = l, size = (800, 300), xlabel = "", ylabel = "" )
xlims!(500, 673)

End of the notebook.