ScoreDrivenModels.jl: a Julia Package for Generalized Autoregressive Score Models
Guilherme Bodin, Raphael Saavedra, Cristiano Fernandes, Alexandre Street
SScoreDrivenModels.jl : a
Julia
Package forGeneralized Autoregressive Score Models
Guilherme Bodin
PUC-RioRio de Janeiro, Brazil
Raphael Saavedra
Invenia LabsCambridge, UK
Cristiano Fernandes
PUC-RioRio de Janeiro, Brazil
Alexandre Street
PUC-RioRio de Janeiro, Brazil
Abstract
Score-driven models, also known as generalized autoregressive score (GAS) models,represent a class of observation-driven time series models. They possess powerful proper-ties, such as the ability to model different conditional distributions and to consider time-varying parameters within a flexible framework. In this paper, we present
ScoreDriven-Models.jl , an open-source
Julia package for modeling, forecasting, and simulating timeseries using the framework of score-driven models. The package is flexible with respect tomodel definition, allowing the user to specify the lag structure and which parameters aretime-varying or constant. It is also possible to consider several distributions, includingBeta, Exponential, Gamma, Lognormal, Normal, Poisson, Student’s t, and Weibull. Theprovided interface is flexible, allowing interested users to implement any desired distribu-tion and parametrization.
Keywords : score-driven models, generalized autoregressive score models, time series models,time-varying parameters, non-Gaussian models.
1. Introduction
Time series models with time-varying parameters have become increasingly popular over theyears due to their advantages in capturing dynamics of series of interest. According to Cox et al. (1981), the mechanism driving parameter dynamics in this general class of models canbe of two types: parameter-driven, as in state space models (Durbin and Koopman 2012;Koopman et al. et al. et al. a r X i v : . [ s t a t . C O ] A ug ScoreDrivenModels.jl : Generalized Autoregressive Score Models in
Julia
Two of the main advantages of the GAS framework are its ability to consider different non-Gaussian distributions and its flexibility with respect to the updating mechanism, which isdetermined by the chosen distribution determines the model updating mechanism. Theseproperties have led GAS models to be applied in numerous fields, such as finance (Har-vey and Thiele 2016; Ayala and Blazsek 2018), actuaries (Neves et al. et al. et al. et al. et al. . This wide range of ap-plications has motivated the development of software packages for this class of models. Forinstance, there are open-source packages in
Python (Taylor 2016), R (Ardia et al. https://timeserieslab.com ), a free software developed by some of the authors of thetheory developed in (Creal et al. Julia (Bezan-son et al.
ScoreDrivenModels.jl (Bodin et al.
Julia ’s main ad-vantages is to avoid the so-called two-language problem, i.e., the dependence on subroutinesimplemented in lower-level languages such as C / C++ or Fortran . Julia achieves this by pro-viding a high-level programming syntax that allows for rapid prototyping and developmentwithout sacrificing computational performance. Thus, by providing an open-source packagecompletely written in
Julia , we facilitate development and contributions by users while alsomaintaining a high level of code transparency. The package allows users to specify a widevariety of GAS models by choosing the conditional distribution, the autoregressive structure,and which parameters are time-varying. Finally, initialization procedures are implemented toturn the estimation process more robust for the case of seasonal time series.The remainder of this paper is organised as follows. Section 2 provides a brief overview of theGAS framework. In Section 3, the
ScoreDrivenModels.jl package is presented, including themodel specification, estimation, forecasting, and simulation. Section 4 presents examples ofapplications to illustrate the use of the package. Conclusions are drawn in Section 5. Finally,the Appendix provides the derivation of the score for each implemented distribution.
2. Score-Driven Models
Let y t ∈ Y ⊆ R denote the dependent variable of interest, f t ∈ P ⊂ R k be a vector oftime-varying parameters, and Y t = { y , . . . , y t } and F t = { f , f , . . . , f t } denote the sets ofavailable information until time t . We assume that y t is generated by the probability densityfunction conditioned on the available information (past data and time-varying parameters)and on the hyperparameter vector θ , which contains the constant parameters. It follows thatthe predictive distribution of y t has a closed form, represented as: p ( y t | F t , Y t − ; θ ) (1) uilherme Bodin, Raphael Saavedra, Cristiano Fernandes, Alexandre Street f t is givenby the following equation, referred to as a GAS( p , q ) mechanism: f t +1 = ω + p X i =1 A i s t − i +1 + q X j =1 B j f t − j +1 , (2)where ω is a vector of constants, coefficient matrices A i and B j have appropriate dimensionsfor i = 1 , . . . , p and j = 1 , . . . , q , and s t = s t ( y t , f t , F t − , Y t − ; θ ) is an appropriate function ofpast data. The unknown coefficients in Eq. (2) are functions of the vector of hyperparameters θ ; that is, ω = ω ( θ ), A i = A i ( θ ), and B j = B j ( θ ). At instant t , the update of the time-varying f t for the next period t + 1 is conducted through Eq. (2), with s t = I − dt | t − · ∇ t , ∇ t = ∂ ln p ( y t | f t , F t − , Y t − ; θ ) ∂f t , (3)where ∇ t is the called the score and I − dt | t − is the scaled Fisher information of the probabilitydensity p ( y t | f t , F t − ; θ ). The scaling coefficient d commonly takes values in { , , } . It isworth mentioning that in the case where d = , it follows that I − t | t − results from the Choleskydecomposition of I − t | t − .As a consequence of the time-varying mechanism for the distribution parameters presented inEq. (2), the conditional distribution of a GAS model is capable of continuously changing basedon the considered data. For instance, if the time series contains occasional volatility spikes,the model can capture this behavior through the time-varying nature of the parameters. Thisproperty is illustrated in Fig. 1.Figure 1: The GAS framework allows the conditional distribution to continuously changebased on the data. In the GAS updating mechanism (2), the parameter f t ∈ P ⊂ R k is sometimes bounded – forexample, in a Normal GAS model where the variance is time-varying, we would have f t = σ t which can only assume positive values by definition. However, in some cases the recursion canlead to updates of f t / ∈ P . A solution is to reparametrize the equations in order to guarantee ScoreDrivenModels.jl : Generalized Autoregressive Score Models in
Julia f t ∈ P for every update. To that end, we follow the procedure described in Creal et al. (2013)and present an example to illustrate it.Let f t be the vector of time-varying parameters of a Normal distribution. From the propertiesof the distribution, it follows that µ t ∈ R and σ t ∈ R + . Let us define a new time-varyingparameter ˜ f t ∈ R k and a map h : P → R k , which we denote the link function in the package.In the case of the Normal distribution, a useful approach is to have an IdentityLink for µ t and a LogLink for σ t as follows: f t = " µ t σ t , ˜ f t = " µ t ln σ t (4)Note that the use of this parametrization will affect the recursion in Eq. (2) as well as thefinal expressions of s t . Thus, let us derive the new recursion for the (2), but this time withthe guarantee that every update of f t respects f t ∈ P . To that end, we also define the inversemap h − ( · ) denoted in the software as the unlink function: f t = h − ( ˜ f t ). Then, we candefine the GAS updating recursion utilizing the parametrization:˜ f t +1 = ω + p X i =1 A i ˜ s t − i +1 + q X j =1 B j ˜ f t − j +1 (5)The unknown coefficients in (5) remain as functions of θ ; that is, ω = ω ( θ ), A i = A i ( θ ), and B j = B j ( θ ). However, this time, we update the linked version of the time-varying parameters˜ f t using (5). It is important to note that the expressions of ˜ s t are different for each scaling d ∈ { , , } . To compute these, we must use the derivative ˙ h of the map h , which is simplyits Jacobian. In our example, h is defined as h " µ t σ t = " h ( µ t ) h ( σ t ) = " ˜ µ t ˜ σ t = " µ t ln σ t (6)and its Jacobian is defined as˙ h " µ t σ t = ∂h ( µ t ) ∂µ t ∂h ( µ t ) ∂σ t ∂h ( σ t ) ∂µ t ∂h ( σ t ) ∂σ t = " σ t (7)Note that ˙ h is always a diagonal matrix. The reparametrized score derivations for differ-ent scalings and different types of maps are presented in Appendix A. Given the followingdefinitions ˙ h ( f t ) = ∂h ( f t ) ∂ ˜ f t (cid:12)(cid:12)(cid:12)(cid:12) f t , J t | t − J > t | t − = I − t | t − , ˜ ∇ t = (cid:16) ˙ h (cid:17) − ∇ t , (8)then the linked scaled score ˜ s t can be computed as follows:˜ s t = (cid:16) ˙ h (cid:17) − ∇ t , for d = 0 , (9)˜ s t = J t | t − ∇ t , for d = 12 , (10)˜ s t = ˙ h I − t | t − ∇ t , for d = 1 . (11) uilherme Bodin, Raphael Saavedra, Cristiano Fernandes, Alexandre Street The vector of hyperparameters θ can be estimated via maximum likelihood:ˆ θ = arg max θ N X t =1 ln p ( y t | f t , F t − ; θ ) (12)Evaluating the log-likelihood function of the GAS model is particularly simple. Given valuesfor the constant parameters θ , the GAS updating equation (2) outputs the conditional dis-tribution at each time period, which generally has a closed form. Thus, it suffices to look atln p ( y t | f t , F t − ; θ ) for a particular value of θ .Evaluating the analytical derivatives needed to obtain the maximum likelihood is a demandingand sometimes impossible task. As a consequence, a common practice is to numerically eval-uate derivatives using global optimization methods such as L-BFGS (Liu and Nocedal 1989)and Nelder-Mead (Nelder and Mead 1965). Depending on the case, constrained optimizationcan also be applied, using, for instance, a Newton interior points method. Forecasting and simulation of future scenarios are among the main goals in time series analysis.In Blasques et al. (2016), details of the procedure for out-of-sample confidence intervals for thetime-varying parameters are discussed. The procedure discussed in Section 4.1 of Blasques et al. (2016) is currently implemented in
ScoreDrivenModels.jl as follows:1. Given ˆ θ T and the filtered state ˆ f T +1 , draw S values y T +1 , . . . , y ST +1 from the estimatedconditional density at T + 1: y T +1 ∼ p ( y T +1 | ˆ f T +1 , ˆ θ T ) for s = 1 , . . . , S .2. Use y T +1 , . . . , y ST +1 and the recursion (2) to obtain the filtered values ˆ f T +2 , . . . , ˆ f ST +2 .3. Repeat steps 1 and 2 H times for H steps ahead generating one new value of y and f per scenario s .Once the procedure is over, S scenarios for the observations within the entire horizon, y sT + k for k = 1 , . . . , H and s = 1 , . . . , S have been simulated. Based on these set of scenario, one cancalculate quantile forecasts, build empirical distributions, or use them to feed decision underuncertainty models, such as stochastic programming. Note that this method solely considersthe uncertainty of innovations. The consideration of uncertainty on both innovations andparameters, as discussed in Blasques et al. (2016), is considered future work for the package.
3. The ScoreDrivenModels.jl Package
ScoreDrivenModels.jl enables users to create and estimate score-driven models and to per-form forecasting and simulation while working purely in
Julia . Its API allows users to choosebetween different distributions, scaling values, lag structures, and optimization methods. Thebasic code structure allows contributors to add new distributions and optimization methods;technical details about adding new features are available in the package documentation. In-stallation of the package is easily conducted using the Julia Package manager:
ScoreDrivenModels.jl : Generalized Autoregressive Score Models in
Julia pkg> add ScoreDrivenModels
To create a
Model , the user must specify 1) the desired distribution, 2) the scaling, 3) the lagstructure, and 4) which parameters should be considered time-varying.1. The lag structure in a GAS( p, q ) model can be specified in two ways: either throughintegers p and q , which results in all lags from to p and to q being added, or througharrays of integers ps and qs containing only the desired lags.2. To specify the distribution, the user needs to choose a distribution among the availableones that have an interface with Distributions.jl . The list of available distributions isdisplayed in Table 1. Furthermore, we refer the interested reader to Appendix B, wherewe provide details on the score calculations for each probability density made availablein the package.3. The scaling is specified by defining the value of d , which can be 0, 1, or , respectivelythe identity scaling, inverse scaling, and inverse square-root scaling.4. In order to define which distribution parameters should be time-varying, the keywordargument time_varying_params can be used. Note that the default behavior is to haveall parameters as time-varying.Distribution Number ofparameters Identityscaling Inversescaling Inversesquare-rootscalingBeta 2 (cid:88) (cid:88) (cid:88) BetaLocationScale 4 (cid:88) – –Exponential 1 (cid:88) (cid:88) (cid:88)
Gamma 2 (cid:88) (cid:88) (cid:88)
LogitNormal 2 (cid:88) (cid:88) (cid:88)
LogNormal 2 (cid:88) (cid:88) (cid:88)
NegativeBinomial 2 (cid:88) – –Normal 2 (cid:88) (cid:88) (cid:88)
Poisson 1 (cid:88) (cid:88) (cid:88)
TDist 1 (cid:88) (cid:88) (cid:88)
TDistLocationScale 3 (cid:88) (cid:88) (cid:88)
Weibull 2 (cid:88) – –Table 1: List of currently implemented distributions and scalings.
NaN within the
Model structure. As an example, a GAS(1 ,
2) model withlognormal distribution and inverse square-root scaling can be created by writing the followingline of code: uilherme Bodin, Raphael Saavedra, Cristiano Fernandes, Alexandre Street julia> Model(1, 2, LogNormal, 0.5) Model{LogNormal,Float64}([NaN, NaN], Dict(1=>[NaN 0.0; 0.0 NaN]),Dict(2=>[NaN 0.0; 0.0 NaN],1=>[NaN 0.0; 0.0 NaN]), 0.5)Dict is the
Julia data structure for dictionaries. Its use allows code flexibility enabling compu-tational simplifications for complex lag structures. As displayed above, the unknown constantparameters to be estimated are set as
NaN . In this case, the constant parameters consideredin vector ω are A , B , and B .In some applications, however, the user might define only one of the distribution parametersas time-varying. In the example below, the only time-varying parameter is µ t , so the keywordargument time_varying_params indicates a vector with only one element, [1] , representingthe first parameter of the lognormal distribution. A table that indicates the distributionparameters and their orders is available in the package documentation. The choice of thetime-varying parameter can be expressed by the following code: julia> Model(1, 2, LogNormal, 0.5; time_varying_params = [1]) Model{LogNormal,Float64}([NaN, NaN], Dict(1=>[NaN 0.0; 0.0 0.0]),Dict(2=>[NaN 0.0; 0.0 0.0],1=>[NaN 0.0; 0.0 0.0]), 0.5)
Users can also specify the lag structure by passing only the lags of interest. Note that thisfeature is equivalent to defining that matrices A i and B j are equal to zero for certain values i and j . An example is a model that uses lags 1 and 12, which means that only the matrices A , A , B , and B have nonzero entries: julia> Model([1, 12], [1, 12], LogNormal, 0.5) Model{Normal,Float64}([NaN, NaN],Dict(12=>[NaN 0.0; 0.0 NaN],1=>[NaN 0.0; 0.0 NaN]),Dict(12=>[NaN 0.0; 0.0 NaN],1=>[NaN 0.0; 0.0 NaN]), 0.5)
Once the model is specified, the next step is estimation. Users can choose from different opti-mization methods provided by
Optim.jl (Mogensen and Riseth 2018). Since this optimizationproblem is non-convex, there is no guarantee that the optimal value found by the optimizationmethod is the global optimum. To increase the chances of finding the global optimum, werun the optimization algorithm for different initial parameter values. The default method isNelder-Mead with 3 random initial parameter values, but the optimization interface is highlyflexible. Users can customize convergence tolerances, choose initial parameter values, and,depending on the optimization method, choose bounds for the parameters. By default, theseinitial values are the unconditional mean of f t +1 which is given by E [ f t +1 ] = ω (cid:18) I − q X j =1 B j (cid:19) − . (13) ScoreDrivenModels.jl : Generalized Autoregressive Score Models in
Julia
As an illustration, let us estimate a GAS model using the same data and specification usedin the R package
GAS (Ardia et al. fit! , the data is alsoavailable in package repository (Bodin et al. julia> Random.seed!(123)julia> y = vec(readdlm("../test/data/cpichg.csv"))julia> gas = Model(1, 1, TDistLocationScale, 0.0,time_varying_params=[1, 2])julia> f = fit!(gas, y)
Round 1 of 3 - Log-likelihood: -178.2064944794775Round 2 of 3 - Log-likelihood: -178.20649327545632Round 3 of 3 - Log-likelihood: -178.20649537930336
Users also have the option to check more detailed results of the optimization procedure bychanging the keyword argument verbose . The default value of this argument is 1; to checkthe optimization summary, users should set the verbose level to 2, and to see the value of theobjective function at each iteration of the optimization, it should be set to 3. To avoid theprinting of outputs, users can set verbose = 0 . To illustrate, results with level 2 is depictedbelow. julia> gas = Model(1, 1, TDistLocationScale, 0.0,time_varying_params=[1, 2])julia> f = fit!(gas, y; verbose=2)
Round 1 of 3 - Log-likelihood: -178.20649402602353Round 2 of 3 - Log-likelihood: -178.20649463073832Round 3 of 3 - Log-likelihood: -178.2064932319776Best initial_point optimization result:* Status: success* Candidate solutionMinimizer: [3.74e-02, -2.60e-01, 1.88e+00, ...]Minimum: 1.782065e+02* Found withAlgorithm: Nelder-MeadInitial Point: [3.13e-02, 1.61e-01, 6.53e-02, ...]* Convergence measuresstandard-deviation <= 1.0e-06* Work counters uilherme Bodin, Raphael Saavedra, Cristiano Fernandes, Alexandre Street Seconds run: 0 (vs limit Inf)Iterations: 1345f(x) calls: 2156
As mentioned before, while the maximization of the log-likelihood is done by default throughthe Nelder-Mead method with 3 random initial values, these features can be changed by theuser. For example, to use the L-BFGS algorithm with 5 random initial values: julia> gas = Model(1, 1, TDistLocationScale, 0.0,time_varying_params=[1, 2])julia> f = fit!(gas, y; opt_method=LBFGS(gas, 5))
Once the estimation step is finished, the user can query the results by calling the function fit_stats : julia> fit_stats(f) --------------------------------------------------------Distribution: Distributions.LocationScale{Float64,TDist{Float64}}Number of observations: 276Number of unknown parameters: 7Log-likelihood: -178.2065AIC: 370.4130BIC: 395.7558--------------------------------------------------------Parameter Estimate Std.Error t stat p-valueomega_1 0.0374 0.0311 1.2016 0.2686omega_2 -0.2599 0.1409 -1.8454 0.1075omega_3 1.8758 0.2914 6.4380 0.0004A_1_11 0.0717 0.0184 3.8884 0.0060A_1_22 0.4538 0.2139 2.1216 0.0715B_1_11 0.9432 0.0272 34.6438 0.0000B_1_22 0.8556 0.0743 11.5141 0.0000 This result matches the example discussed in Ardia et al. (2019) with the exception of omega_3 , due to a difference in parametrization between the two packages. Once the pa-rameter is recovered to its original parametrization, the result becomes the same.
Forecasting in this framework is done by simulation as per Section 2.4. Function forecast runs the procedure proposed by Blasques et al. (2016) and returns a
Forecast structure thathas the expected value for time-varying parameters, observations, and the related scenariosused to find them. By default the structure also stores the 2.5%, 50% and 97.5% quantiles.Next, we will present forecasting results using the previously estimated US inflation data. Inthe example below, the first column is the location parameter, the second column is the scaleparameter, and the third column represents the degrees of freedom parameter.0
ScoreDrivenModels.jl : Generalized Autoregressive Score Models in
Julia julia> forec = forecast(y, gas, 12)julia> forec.parameter_forecast
We can also obtain the scenarios of observations, y sT + k for k = 1 , . . . , H and s = 1 , . . . , S ,that generated the above forecasted values as follows: julia> forec.observation_scenarios For the sake of clarity, the forecast ˆ y T + k is the sample average of the scenarios of observations,i.e., ˆ y T + k = S − P Ss =1 y sT + k .
4. Applications
The Brazilian system operator regularly publishes an indicator of how much energy can begenerated from water inflows among all hydropower plants in the country. This indicator is uilherme Bodin, Raphael Saavedra, Cristiano Fernandes, Alexandre Street µ t . Finally, ANE scenarios are simulated based on the fitted hyperparameters.Figure 2: ANE in the Northeastern region of BrazilIn this model we utilize lags 1, 2, 11 and 12 for the score and autoregressive components; lags1 and 2 capture the short-term trend of the series, while lags 11 and 12 capture the seasonalaspect of water inflow. The model can be written as follows: µ t = ω + A s t − + A s t − + A s t − + A s t − + B µ t − + B µ t − + B µ t − + B µ t − ln( σ t ) = ω (14)In this model we have also considered that the initial parameters µ , . . . µ , σ , . . . , σ , s , . . . s are calculated heuristically through maximum likelihood estimation in each sea-sonal component of the model. This way, the parameters µ , σ are calculated by fitting a2 ScoreDrivenModels.jl : Generalized Autoregressive Score Models in
Julia lognormal distribution in the observations y , y , y , . . . , y n . This procedure was used inHoeltgebaum et al. (2018). Note that this procedure, which is implemented by the function dynamic_initial_params in our package, is relevant for ensuring good estimation resultswhen considering seasonal time series. The model is estimated as follows: julia> Random.seed!(123)julia> y = vec(readdlm("../test/data/ane_northeastern.csv"))julia> y_train = y[1:400]julia> gas = Model([1, 2, 11, 12], [1, 2, 11, 12], LogNormal, 0.0;+ time_varying_params=[1]julia> initial_params = dynamic_initial_params(y_train, gas)julia> f = ScoreDrivenModels.fit!(gas, y_train; initial_params=initial_params)julia> fit_stats(f) --------------------------------------------------------Distribution: LogNormalNumber of observations: 400Number of unknown parameters: 10Log-likelihood: -779.7883AIC: 1579.5766BIC: 1619.4912--------------------------------------------------------Parameter Estimate Std.Error t stat p-valueomega_1 0.0135 0.0282 0.4807 0.6411omega_2 -2.8408 0.0720 -39.4651 0.0000A_1_11 -0.0378 0.0051 -7.4109 0.0000A_2_11 0.0047 0.0034 1.4052 0.1903A_11_11 -0.0178 0.0045 -3.9838 0.0026A_12_11 0.0576 0.0052 11.0468 0.0000B_1_11 -0.4784 0.0444 -10.7815 0.0000B_2_11 0.4682 0.0449 10.4210 0.0000B_11_11 -0.3055 0.0928 -3.2933 0.0081B_12_11 1.3088 0.0962 13.6080 0.0000 Once the model is estimated we can generate 1000 ANE scenarios for the next 60 monthsfollowing the methodology discussed in Section 2.4. We will utilize forecast in order toobtain the quantiles as well. The data set used to estimate the model used data from January1961 until April 1994. To illustrate the adequacy of our model forecast, we present in Fig. 3an out-of-sample study, where the simulated scenarios (gray lines), and related quantiles (reddotted lines) from May 1994 until April 1999 are contrasted to the actual observed values. julia> forec = forecast(y_train, gas, 60;S = 1_000, initial_params = initial_params)
One of the advanced features of
ScoreDrivenModels.jl is allowing users to change the defaultparametrization. An example of a different parametrization is the GARCH(1, 1) GAS model. uilherme Bodin, Raphael Saavedra, Cristiano Fernandes, Alexandre Street h is theidentity (we refer the interested reader to Appendix A.3 for further details). To ensure theequivalence, the user must define the identity link function for both time-varying parameters.To do this the user must overwrite three ScoreDrivenModel.jl methods as detailed in thefollowing example: julia> function ScoreDrivenModels.link!(param_tilde::Matrix{T},::Type{Normal}, param::Matrix{T}, t::Int) where Tparam_tilde[t, 1] = link(IdentityLink, param[t, 1])param_tilde[t, 2] = link(IdentityLink, param[t, 2])returnendjulia> function ScoreDrivenModels.unlink!(param::Matrix{T}, ::Type{Normal},param_tilde::Matrix{T}, t::Int) where Tparam[t, 1] = unlink(IdentityLink, param_tilde[t, 1])param[t, 2] = unlink(IdentityLink, param_tilde[t, 2])returnendjulia> function ScoreDrivenModels.jacobian_link!(aux::AuxiliaryLinAlg{T},::Type{Normal}, param::Matrix{T}, t::Int) where Taux.jac[1] = jacobian_link(IdentityLink, param[t, 1]) ScoreDrivenModels.jl : Generalized Autoregressive Score Models in
Julia aux.jac[2] = jacobian_link(IdentityLink, param[t, 2])returnend
Once the methods have been overwritten for the Normal distribution, the recursion willapply the
IdentityLink in both parameters and the user can proceed to the estimationstep. Below, we run an example provided in Broda and Paolella (2020) for daily Germanmark/British pound exchange rates. Note that there are also bounds being provided for thehyperparameter estimation through lb and ub . julia> y = vec(readdlm("../test/data/BG96.csv"))julia> initial_params = [mean(y) var(y)]julia> ub = [1.0, 1.0, 0.5, 1.0]julia> lb = [-1.0, 0.0, 0.0, 0.5]julia> gas = Model(1, 1, Normal, 1.0, time_varying_params = [2])julia> initial_point = [0.0, 0.5, 0.25, 0.75]julia> f = fit!(gas, y; initial_params = initial_params,opt_method = IPNewton(gas, [initial_point]ub=ub, lb=lb))Round 1 of 1 - Log-likelihood: -1106.598367006442julia> fit_stats(f) --------------------------------------------------------Distribution: NormalNumber of observations: 1974Number of unknown parameters: 4Log-likelihood: -1106.5984AIC: 2221.1967BIC: 2243.5480--------------------------------------------------------Parameter Estimate Std.Error t stat p-valueomega_1 -0.0062 0.0085 -0.7373 0.5019omega_2 0.0108 0.0029 3.7732 0.0196A_1_22 0.1534 0.0266 5.7726 0.0045B_1_22 0.9593 0.0144 66.6015 0.0000 The results obtained above are the same as the ones found in the Section 4.1 of Broda andPaolella (2020), with the exception of one parameter. This is due to a difference in theparametrizations of the models – the demonstration of the equivalence between the hyper-parameters of the GAS model and the GARCH model can be found in Appendix A.3. Notethat, following the demonstration in A.3, we have β = B − A = 0 . rugarch (Ghalanos 2020) and ob-tained the same result, thus further illustrating the equivalence of the GARCH(1, 1) and theNormal GAS(1, 1) under this parametrization. uilherme Bodin, Raphael Saavedra, Cristiano Fernandes, Alexandre Street
5. Conclusion
The
ScoreDrivenModels.jl package provides a general framework for score-driven models andrepresents a user-friendly, off-the-shelf tool. The package is fully implemented in Julia, thusnot depending on subroutines written in lower-level languages such as C or Fortran . Themodel specification is flexible and allows defining any desired lag structure, as well as anydistribution due to a tailored dependency on the
Distributions.jl package. The estimationprocedure is based on numerical optimization algorithms such as Nelder-Mead or L-BFGS andemploys the well-known package
Optim.jl . Special initialization procedures are implementedto robustify the estimation process for the case of seasonal time series. Available forecastingand simulation procedures allow users to study future data from the estimated model. Finally,the examples provided in Section 4 illustrate the functionalities of the package as well aspossible applications. The software continues to evolve, new features such as an heuristicfor the initial hyperparameter and unobserved components modeling are considered as futureresearch topics. The software documentation can be found in https://lampspuc.github.io/ScoreDrivenModels.jl/latest/
Acknowledgements
This study was funded in part by the Coordenação de Aperfeiçoamento de Pessoal de NívelSuperior - Brasil (CAPES) - Finance Code 001, by the Conselho Nacional de Desenvolvi-mento Científico e Tecnológico (CNPq), and by the Energisa Group through the R&D projectANEEL PD-00405-1701/2017.
References
Ardia D, Boudt K, Catania L (2019). “Generalized autoregressive score models in R: TheGAS package.”
Journal of Statistical Software , (1). ISSN 15487660. doi:10.18637/jss.v088.i06 . .Ayala A, Blazsek S (2018). “Score-driven copula models for portfolios of two risky assets.” TheEuropean Journal of Finance , (18), 1861–1884. doi:10.1080/1351847X.2018.1464488 .Bezanson J, Edelman A, Karpinski S, Shah V (2017). “Julia: A fresh approach to numericalcomputing.” SIAM review , (1), 65–98. doi:10.1137/141000671 .Blasques F, Koopman SJ, Łasak K, Lucas A (2016). “In-sample confidence bands and out-of-sample forecast bands for time-varying parameters in observation-driven models.” Interna-tional Journal of Forecasting , (3), 875–887. doi:10.1016/j.ijforecast.2015.11.018 .Bodin G, Saavedra R, Fernandes C (2020). “ ScoreDrivenModels.jl [Online]. Available: https://github.com/LAMPSPUC/ScoreDrivenModels.jl .”Bollerslev T (1986). “Generalized autoregressive conditional heteroskedasticity.”
Journal ofeconometrics , (3), 307–327.Broda SA, Paolella MS (2020). “ARCHModels.jl: Estimating ARCH models in Julia.”6 ScoreDrivenModels.jl : Generalized Autoregressive Score Models in
Julia
Cox DR, Gudmundsson G, Lindgren G, Bondesson L, Harsaae E, Laake P, Juselius K, Lau-ritzen SL (1981). “Statistical Analysis of Time Series: Some Recent Developments [withDiscussion and Reply].”
Scandinavian Journal of Statistics , (2), 93–115. ISSN 03036898,14679469. URL .Creal D, Koopman SJ, Lucas A (2013). “Generalized autoregressive score models withapplications.” Journal of Applied Econometrics , (5), 777–795. ISSN 08837252. doi:10.1002/jae.1279 .de Melo MAB, Fernandes CA, de Melo EF (2018). “Forecasting aggregate claims using score-driven time series models.” Statistica Neerlandica , (3), 354–374. doi:10.1111/stan.12139 .Durbin J, Koopman SJ (2012). Time Series Analysis by State Space Methods . Oxford Uni-versity Press. doi:10.1093/acprof:oso/9780199641178.001.0001 .Engle RF, Russell JR (1998). “Autoregressive Conditional Duration: A New Model forIrregularly Spaced Transaction Data.”
Econometrica , (5), 1127–1162. ISSN 00129682,14680262. URL .Ghalanos A (2020). “rugarch: Univariate GARCH Models. R package version 1.4-2 [Online].Available: https://CRAN.R-project.org/package=rugarch .”Harvey AC (2013). Dynamic Models for Volatility and Heavy Tails: With Applications toFinancial and Economic Time Series . Econometric Society Monographs. Cambridge Uni-versity Press. doi:10.1017/CBO9781139540933 .Harvey AC, Thiele S (2016). “Testing against changing correlation.”
Journal of EmpiricalFinance , (Part B), 575–589. ISSN 09275398. doi:10.1016/j.jempfin.2015.09.003 .Hoeltgebaum H, Fernandes C, Street A (2018). “Generating joint scenarios for renewablegeneration: The case for non-gaussian models with time-varying parameters.” IEEE Trans-actions on Power Systems , (6), 7011–7019. ISSN 08858950. doi:10.1109/TPWRS.2018.2838050 .Koopman SJ, Harvey AC, Doornik JA, Shephard N (2000). “STAMP 6.0: Structural TimeSeries Analyser, Modeller and Predictor.” London: Timberlake Consultants .Liu DC, Nocedal J (1989). “On the limited memory BFGS method for large scale optimiza-tion.”
Mathematical programming , (1-3), 503–528.Mogensen PK, Riseth AN (2018). “Optim: A mathematical optimization package for Julia.” Journal of Open Source Software , (24). doi:10.21105/joss.00615 .Nani A, Gamoudi I, El Ghourabi M (2019). “Value-at-risk estimation by LS-SVR and FS-LS-SVR based on GAS model.” Journal of Applied Statistics , (12), 2237–2253. doi:10.1080/02664763.2019.1584161 .Nelder JA, Mead R (1965). “A Simplex Method for Function Minimization.” The Com-puter Journal , (4), 308–313. ISSN 0010-4620. doi:10.1093/comjnl/7.4.308 . https://academic.oup.com/comjnl/article-pdf/7/4/308/1013182/7-4-308.pdf , URL https://doi.org/10.1093/comjnl/7.4.308 . uilherme Bodin, Raphael Saavedra, Cristiano Fernandes, Alexandre Street Insurance: Mathematicsand Economics , , 48–57. ISSN 01676687. doi:10.1016/j.insmatheco.2017.04.004 .Patton AJ, Ziegel JF, Chen R (2019). “Dynamic semiparametric models for expected shortfall(and value-at-risk).” Journal of Econometrics , (2), 388–413.Saavedra R (2017). “A study on the impact of El Niño Southern Oscillation on hydro powergeneration in Brazil.” doi:10.17771/PUCRio.acad.32290 .Saavedra R, Bodin G, Souto M (2019). “StateSpaceModels.jl: a Julia Package for Time-SeriesAnalysis in a State-Space Framework.” arXiv preprint arXiv:1908.01757 .Taylor R (2016). “PyFlux.” [Online]. Available: https://github.com/RJT1990/pyflux . A. Parametrizations
The use of different link functions and scaling coefficient give rise to different probabilityfunction parametrizations and models. Thus, in this appendix we explore some modelingvariants due to the choice of link functions and values for the scaling coefficient.
A.1. Possible Parametrizations
There are three mapping functions available in
ScoreDrivenModels.jl . Within the contextof the software we call them
Links . Each one of them can be used as the default mappingfunction for a given distribution implemented in the package.• Identity link: ˜ f = f where f ∈ R and ˜ f ∈ R • Log link: ˜ f = ln ( f − a ) where f ∈ [ a, ∞ ) and ˜ f ∈ R • Logit link: ˜ f = ln (cid:16) f − ab − f (cid:17) where f ∈ [ a, b ] and ˜ f ∈ R A.2. Score Derivations for Different Scaling Values
In this subsection we derive expressions (9), (10), and (11).
Scaling d = 0
For d = 0 the score is simply equal to˜ ∇ t = ∂ ln p ( y t | y t − , f t ) ∂ ˜ f t , (15)which is equivalent to ˜ ∇ t = ∂f t ∂ ˜ f t · ∂ ln p ( y t | y t − , f t ) ∂f t . (16)Notice that one can show by the inverse function theorem that this is the inverse of theJacobian ˙ h . Thus, it follows that ˜ ∇ t = (cid:16) ˙ h (cid:17) − ∇ t . (17)8 ScoreDrivenModels.jl : Generalized Autoregressive Score Models in
Julia
By proceeding the calculus above indicated, we show that˜ ∇ t = ∂ ln p ( y t | y t − ,µ,σ ) ∂µ∂ ln p ( y t | y t − ,µ,σ ) ∂ ˜ σ = " σ t − ∂ ln p ( y t | y t − ,µ,σ ) ∂µ∂ ln p ( y t | y t − ,µ,σ ) ∂σ (18)˜ ∇ t = ∂ ln p ( y t | y t − ,µ,σ ) ∂ ˜ µ∂ ln p ( y t | y t − ,µ,σ ) ∂ ˜ σ = " σ t ∂ ln p ( y t | y t − ,µ,σ ) ∂µ∂ ln p ( y t | y t − ,µ,σ ) ∂σ (19)˜ ∇ t = ∂ ln p ( y t | y t − ,µ,σ ) ∂ ˜ µ∂ ln p ( y t | y t − ,µ,σ ) ∂ ˜ σ = ∂ ln p ( y t | y t − ,µ,σ ) ∂µ∂ ln p ( y t | y t − ,µ,σ ) ∂σ · σ t . (20) Scaling d = 1/2
The original scaled score is s t = J t | t − ∇ t as in considered in Creal et al. (2013). Thus, recallthat I − t | t − = J t | t − J > t | t − . (21)In this case, the new scaled score is ˜ s t = ˜ J t | t − ˜ ∇ t . Thus, if ˜ ∇ t = (cid:16) ˙ h (cid:17) − ∇ t and ˜ J t | t − is yetto be calculated, we have I t | t − = E h ∇ t ∇ > t i (22)˜ I t | t − = E (cid:20)(cid:16) ˙ h (cid:17) − ∇ t ∇ > t (cid:16) ˙ h (cid:17) − (cid:21) (23)˜ I t | t − = (cid:16) ˙ h (cid:17) − E h ∇ t ∇ > t i (cid:16) ˙ h (cid:17) − (24)˜ I t | t − = (cid:16) ˙ h (cid:17) − I t | t − (cid:16) ˙ h (cid:17) − . (25)As (cid:16) ˙ h (cid:17) − is a diagonal matrix, in the above development we omitted the transpose operator.So, the inverse of the reparametrized information matrix is equal to˜ I − t | t − = ˙ h I − t | t − ˙ h (26)˜ I − t | t − = ˙ h J t | t − J > t | t − ˙ h. (27)Hence, it follows that ˜ J t | t − ˜ J > t | t − = ˙ h J t | t − J > t | t − ˙ h (28)˜ J t | t − = ˙ h J t | t − , (29)which lead us to conclude that for this type of scaling ˜ s t = s t . This becomes clear in thefollowing development: ˜ s t = ˜ J t | t − ˜ ∇ t (30)˜ s t = ˙ h J t | t − (cid:16) ˙ h (cid:17) − ∇ t , (31)In the case that J t | t − is diagonal, because ˙ h is diagonal we have˜ s t = J t | t − ∇ t = s t . (32) uilherme Bodin, Raphael Saavedra, Cristiano Fernandes, Alexandre Street Scaling d = 1
In this case, the scaled score is equal to s t = I − t | t − ∇ t and the reparametrized scaled score isequal to ˜ s t = ˜ I − t | t − ˜ ∇ t . As previously calculated we have that˜ I − t | t − = ˙ h I − t | t − ˙ h (33)˜ ∇ t = (cid:16) ˙ h (cid:17) − ∇ t . (34)Therefore, ˜ s t = ˜ I − t | t − ˜ ∇ t (35)˜ s t = ˙ h I − t | t − ˙ h (cid:16) ˙ h (cid:17) − ∇ t (36)˜ s t = ˙ h I − t | t − ∇ t (37)˜ s t = ˙ hs t . (38) A.3. Different Parametrizations Lead to Different Models
One of the examples given in Creal et al. (2013) shows the equivalence between a GARCH(1,1) and GAS(1, 1) with Normal distribution and d = 1. An important note on this fact is thatthe models are only equivalent if the variance, σ , is considered a time-varying parameter,i.e., if the link is the identity function. Therefore, if a log link is used to ensure a positivevalue for σ , for instance, the equivalence will not hold. As this calculations are not shown inany other work and they reveal relevant insights that can be tested in through our software,we provide further details in the sequel.Thus, let us develop the GAS recursion for both cases. Recall the Normal probability densityfunction (pdf): p ( y t | y t − , µ, σ ) = 1 √ πσ e (cid:16) − ( yt − µ )22 σ (cid:17) . (39)The log-pdf is: ln p ( y t | y t − , µ, σ ) = −
12 ln 2 π −
12 ln σ −
12 ( y t − µ ) σ . (40)To calculate the score with respect to σ we need to calculate the following derivative: ∂ ln p ( y t | y t − , µ, σ ) ∂σ = − σ + ( y t − µ ) σ . (41)Then, we calculate the Fisher information as follows: − E " ∂ ln p ( y t | y t − , µ, σ ) ∂σ ∂σ = 12 σ . (42)Now if we proceed to write the GAS(1,1) recursion using the inverse scaling d = 1, we find σ t +1 = ω + A s t + B σ t (43) σ t +1 = ω + A (( y t − µ t ) − σ t ) + B σ t . (44)0 ScoreDrivenModels.jl : Generalized Autoregressive Score Models in
Julia
By assuming µ = 0 under correct specification, the recursion becomes σ t +1 = ω + A ( y t − σ t ) + B σ t , (45)which is equivalent to the GARCH(1,1) model σ t +1 = α + α y t + β σ t . (46)Note that there is sufficient degrees of freedom to make the correspondence between theparameters of the two models, e.g., α = ω, α = A , β = B − A .Now let us work with a different parametrization to assure a positive value to σ . Theapproach suggested in Creal et al. (2013) is to use a map h ( · ) = ln( · ), i.e., ˜ σ = ln σ . Whenwe use this parametrization the recursion adapts the following way: ( σ t = h − (˜ σ t ) , ˜ σ t +1 = ω + A ˜ s t + B ˜ σ t . (47)And we shown in the previous section that ˜ s t = ˙ hs t . In this case ˙ h = σ t , so the recursionassumes the following form: ( σ t = h − (˜ σ t ) , ˜ σ t +1 = ω + A y t − σ t σ t + B ˜ σ t . (48)If we rewrite the recursion solely in terms of σ we haveln (cid:16) σ t +1 (cid:17) = ω + A y t − σ t σ t + B ln (cid:16) σ t (cid:17) . (49)In this case however, it is impossible to choose the parameters ω, A and B to meet a recursionequivalent to (46). B. Score Calculations
GAS models can be still considered a relatively recent technology. In this context, this paperalso aims to provide users with a technical reference for the
ScoreDrivenModels.jl software.Therefore, in this Appendix, we provide detailed information on the score calculation for eachdistribution considered in the package.
B.1. Beta
Density function p ( y t | y t − , α, β ) = y α − t (1 − y t ) β − B ( α, β ) where B ( α, β ) = Γ( α )Γ( β )Γ( α + β ) ,α ∈ R + , β ∈ R + , Γ( · ) is the Gamma function. E [ y t | y t − ] = αα + β , VAR [ y t | y t − ] = αβ ( α + β ) ( α + β + 1) uilherme Bodin, Raphael Saavedra, Cristiano Fernandes, Alexandre Street ∇ t = " ∂ ln p ( y t | y t − ,α,β ) ∂α∂ ln p ( y t | y t − ,α,β ) ∂β ln p ( y t | y t − , a, c, α, β ) = ( α −
1) ln y t + ( β −
1) ln(1 − y t ) − ln B ( α, β ) ∇ αt = ∂ ln p ( y t | y t − , α, β ) ∂α = ln y t + ψ ( α + β ) − ψ ( α ) ∇ βt = ∂ ln p ( y t | y t − , a, c, α, β ) ∂β = ln(1 − y t ) + ψ ( α + β ) − ψ ( β ) B.2. Beta location scale
Density function p ( y t | y t − , a, c, α, β ) = ( y t − a ) α − ( c − y t ) β − ( c − a ) B ( α, β ) where B ( α, β ) = Γ( α )Γ( β )Γ( α + β ) ,a, c ∈ R α, β ∈ R + E [ y t | y t − ] = a + ( c − a ) αα + β , VAR [ y t | y t − ] = ( c − a ) αβ ( α + β ) ( α + β + 1)Score calculation ∇ t = ∂ ln p ( y t | y t − ,a,c,α,β ) ∂a∂ ln p ( y t | y t − ,a,c,α,β ) ∂c∂ ln p ( y t | y t − ,a,c,α,β ) ∂α∂ ln p ( y t | y t − ,a,c,α,β ) ∂β ln p ( y t | y t − , a, c, α, β ) = ( α −
1) ln ( y t − a ) + ( β −
1) ln( c − y t ) − ( α + β −
1) ln ( c − a ) − ln B ( α, β ) ∇ at = ∂ ln p ( y t | y t − , a, c, α, β ) ∂a = − α + 1 y − a + α + β − c − a ∇ ct = ∂ ln p ( y t | y t − , a, c, α, β ) ∂c = β − c − y − α + β − c − a ∇ αt = ∂ ln p ( y t | y t − , a, c, α, β ) ∂α = ln( y t − a ) − ln( c − a ) + ψ ( α + β ) − ψ ( α ) ∇ βt = ∂ ln p ( y t | y t − , a, c, α, β ) ∂β = ln( c − y t ) − ln( c − a ) + ψ ( α + β ) − ψ ( β )2 ScoreDrivenModels.jl : Generalized Autoregressive Score Models in
Julia
B.3. Exponential
Density function p ( y t | y t − , λ ) = λe − λy t λ ∈ R + E [ y t | y t − ] = 1 λ , VAR [ y t | y t − ] = 1 λ Score calculation ∇ t = h ∂ ln p ( y t | y t − ,λ ) ∂λ i ln p ( y t | y t − , λ, k ) = ln λ − e − λy t ∇ λt = ∂ ln p ( y t | y t − , λ ) ∂λ = 1 λ − y t B.4. Gamma
Density function p ( y t | y t − , α, k ) = y α − t e − ytk Γ( α ) k α , α ∈ R + , k ∈ R + E [ y t | y t − ] = αk, VAR [ y t | y t − ] = αk Score calculation ∇ t = " ∂ ln p ( y t | y t − ,α,k ) ∂α∂ ln p ( y t | y t − ,α,k ) ∂k ln p ( y t | y t − , α, k ) = ( α −
1) ln y t − y t k − ln Γ( α ) − α ln k ∇ αt = ∂ ln p ( y t | y t − , α, k ) ∂α = ln y t − ψ ( α ) − ln k ∇ kt = ∂ ln p ( y t | y t − , α, k ) ∂k = y t k − αk B.5. Logit-Normal
Density function p ( y t | y t − , µ, σ ) = 1 y t (1 − y t ) √ πσ e (cid:16) − (logit( yt ) − µ )22 σ (cid:17) , µ ∈ R , σ ∈ R + E [ y t | y t − ] = e (cid:16) µ + σ (cid:17) , VAR [ y t | y t − ] = (cid:16) e σ − (cid:17) e ( µ + σ ) uilherme Bodin, Raphael Saavedra, Cristiano Fernandes, Alexandre Street ∇ t = ∂ ln p ( y t | y t − ,µ,σ ) ∂µ∂ ln p ( y t | y t − ,µ,σ ) ∂σ ln p ( y t | y t − , µ, σ ) = − ln ( y t (1 − y t )) −
12 ln 2 πσ − (logit ( y t ) − µ ) σ ∇ µt = ∂ ln p ( y t | y t − , µ, σ ) ∂µ = logit ( y t ) − µσ ∇ σ t = ∂ ln p ( y t | y t − , µ, σ ) ∂σ = − σ − logit ( y t ) − µ ) σ ! B.6. Lognormal
Density function p ( y t | y t − , µ, σ ) = 1 y t √ πσ e (cid:16) − (ln yt − µ )22 σ (cid:17) , µ ∈ R , σ ∈ R + E [ y t | y t − ] = e (cid:16) µ + σ (cid:17) , VAR [ y t | y t − ] = (cid:16) e σ − (cid:17) e ( µ + σ )Score calculation ∇ t = ∂ ln p ( y t | y t − ,µ,σ ) ∂µ∂ ln p ( y t | y t − ,µ,σ ) ∂σ ln p ( y t | y t − , µ, σ ) = − ln y t −
12 ln 2 πσ − (ln y t − µ ) σ ∇ µt = ∂ ln p ( y t | y t − , µ, σ ) ∂µ = ln y t − µσ ∇ σ t = ∂ ln p ( y t | y t − , µ, σ ) ∂σ = − σ − (ln y t − µ ) σ ! B.7. Negative binomial
Density function p ( y t | y t − , r, p ) = Γ( y t + r ) y t !Γ( r ) p r (1 − p ) y t , r ∈ R + , p ∈ [0 , E [ y t | y t − ] = pr − p , VAR [ y t | y t − ] = pr (1 − p ) ScoreDrivenModels.jl : Generalized Autoregressive Score Models in
Julia
Score calculation ∇ t = " ∂ ln p ( y t | y t − ,r,p ) ∂r∂ ln p ( y t | y t − ,r,p ) ∂p ln p ( y t | y t − , r, p ) = ln Γ( y t + r ) − (ln y t !Γ( r )) + r ln p + y t ln(1 − p ) ∇ rt = ∂ ln p ( y t | y t − , r, p ) ∂r = ψ ( y t + r ) − ψ ( r ) + ln p ∇ pt = ∂ ln p ( y t | y t − , r, p ) ∂p = rp − k − p B.8. Normal
Density function p ( y t | y t − , µ, σ ) = 1 √ πσ e (cid:16) − ( yt − µ )22 σ (cid:17) , µ ∈ R , σ ∈ R + E [ y t | y t − ] = µ, VAR [ y t | y t − ] = σ Score calculation ∇ t = ∂ ln p ( y t | y t − ,µ,σ ) ∂µ∂ ln p ( y t | y t − ,µ,σ ) ∂σ ln p ( y t | y t − , µ, σ ) = −
12 ln 2 πσ − ( y t − µ ) σ ∇ µt = ∂ ln p ( y t | y t − , µ, σ ) ∂µ = y t − µσ ∇ σ t = ∂ ln p ( y t | y t − , µ, σ ) ∂σ = − σ − ( y t − µ ) σ ! B.9. Poisson
Density function p ( y t | y t − , λ ) = e − λ λ y t y t ! , λ ∈ R + E [ y t | y t − ] = λ, VAR [ y t | y t − ] = λ Score calculation ∇ t = ∂ ln p ( y t | y t − , λ ) ∂λ ln p ( y t | y t − , λ ) = − λ + y t ln λ − ln y t ! ∇ λt = ∂ ln p ( y t | y t − , λ ) ∂λ = y t − λλ uilherme Bodin, Raphael Saavedra, Cristiano Fernandes, Alexandre Street B.10. Student’s t Density function p ( y t | y t − , ν ) = 1 √ νB (cid:16) , ν (cid:17) y t ν ! − ν +12 , ν ∈ R + E [ y t | y t − ] = 0 , VAR [ y t | y t − ] = σ νν − ν > ∞ for 1 < ν ≤
2, undefined otherwise.Score calculation ∇ t = h ∂ ln p ( y t | y t − ,ν ) ∂ν i ln p ( y t | y t − , ν ) = −
12 ln ν − ln B (cid:18) , ν (cid:19) − (cid:18) ν + 12 (cid:19) ln y t ν ! ∇ νt = ∂ ln p ( y t | y t − ν ) ∂ν = 12 ( ν + 1) y t νy t + ν − ν − ln y t ν + 1 ! + ψ (cid:18) ν + 12 (cid:19) − ψ (cid:18) ν (cid:19)! B.11. Student’s t with Location and Scale Density function p ( y t | y t − , µ, σ , ν ) = 1 √ σ νB (cid:16) , ν (cid:17) y t − µ ) σ ν ! − ν +12 , µ ∈ R , σ ∈ R + ν ∈ R + E [ y t | y t − ] = µ, VAR [ y t | y t − ] = σ νν − ν > ∞ for 1 < ν ≤
2, undefined otherwise.6
ScoreDrivenModels.jl : Generalized Autoregressive Score Models in
Julia
Score calculation ∇ t = ∂ ln p ( y t | y t − ,µ,σ ,ν ) ∂µ∂ ln p ( y t | y t − ,µ,σ ,ν ) ∂σ ∂ ln p ( y t | y t − ,µ,σ ,ν ) ∂ν ln p ( y t | y t − , µ, σ , ν ) = −
12 ln νσ − ln B (cid:18) , ν (cid:19) − (cid:18) ν + 12 (cid:19) ln y t − µ ) σ ν ! ∇ µt = ∂ ln p ( y t | y t − , µ, σ , ν ) ∂µ = ( ν + 1) ( y t − µ )( y t − µ ) + σ ν ∇ σ t = ∂ ln p ( y t | y t − , µ, σ , ν ) ∂σ = − ν σ − ( y t − µ ) σ (cid:16) νσ + ( y t − µ ) (cid:17) ∇ νt = ∂ ln p ( y t | y t − , µ, σ , ν ) ∂ν = 12 ( ν + 1) ( y t − µ ) ν ( y t − µ ) + σ ν − ν − ln ( y t − µ ) σ ν + 1 !! +12 (cid:18) ψ (cid:18) ν + 12 (cid:19) − ψ (cid:18) ν (cid:19)(cid:19) B.12. Weibull
Density function p ( y t | y t − , λ, k ) = kλ (cid:0) y t λ (cid:1) k − e ( − ytλ ) k x ≥ , x < , λ ∈ R + , k ∈ R + E [ y t | y t − ] = λ Γ (1 + 1 /k ) , VAR [ y t | y t − ] = λ " Γ (cid:18) k (cid:19) − (cid:18) Γ (cid:18) k (cid:19)(cid:19) Score calculation ∇ t = " ∂ ln p ( y t | y t − ,λ,k ) ∂k∂ ln p ( y t | y t − ,λ,k ) ∂λ ln p ( y t | y t − , λ, k ) = ln k + ( k −
1) ln y t − k ln λ − (cid:18) y t λ (cid:19) k ∇ λt = ∂ ln p ( y t | y t − , λ, k ) ∂λ = kλ (cid:18) y t λ (cid:19) k − ! ∇ kt = ∂ ln p ( y t | y t − , λ, k ) ∂k = 1 k + ln (cid:18) y t λ (cid:19) − (cid:18) y t λ (cid:19) k ! uilherme Bodin, Raphael Saavedra, Cristiano Fernandes, Alexandre Street C. Computational details
The results in this paper were obtained using
Julia
ScoreDrivenModels v0.1.2.Figure 3 was generated with
Plots.jl v1.4.3. In order to guarantee that the results of theexamples are exactly the same as the ones reported, we recommend the interested user to activate the
Project.toml in the repository’s example folder using Julia’s
Pkg manager.
Affiliation:
Guilherme Bodin (corresponding author)Department of Electrical EngineeringPontifícia Universidade Católica do Rio de JaneiroR. Marquês de São Vicente, 225Gávea, Rio de Janeiro - RJ, BrazilE-mail: [email protected]
Raphael SaavedraInvenia Labs95 Regent StCambridge CB2 1AW, United KingdomE-mail: [email protected]
Cristiano FernandesDepartment of Electrical EngineeringPontifícia Universidade Católica do Rio de JaneiroR. Marquês de São Vicente, 225Gávea, Rio de Janeiro - RJ, BrazilE-mail: [email protected]