Beta Autoregressive Fractionally Integrated Moving Average Models
Guilherme Pumi Marcio Valk Cleber Bisognin Fábio Mariano Bayer Taiane Schaedler Prass

empty
Beta autoregressive fractionally integrated moving average models


 Guilherme Pumii^a,[1]Corresponding author.[-.3cms^aMathematics and Statistics Institute - Universidade Federal do Rio Grande do Sul - 9500,  Bento Gonçalves Avenue - 91509-900, Porto Alegre - RS - Brazil.
], Marcio Valki^a, Cleber Bisognini^b, Fábio Mariano Bayers^b[-.3cms^b Departamento de Estatística and LACESM - Universidade Federal de Santa Maria, Santa Maria - RS - Brazil.
]  and Taiane Schaedler Prasss^a  

[E-mails: guilherme.pumi@ufrgs.br, marciovalk@gmail.com, cbisognin@ufrgs.br, bayer@ufsm.br and taianeprass@gmail.com]

.3cm







In this work we introduce the class of beta autoregressive fractionally integrated moving average models for continuous random variables taking values in the continuous unit interval (0,1).
The proposed model accommodates a set of regressors and a long-range dependent time series structure. We derive the partial likelihood estimator for the parameters of the proposed model, obtain the associated score vector and Fisher information matrix.
We also prove the consistency and asymptotic normality of the estimator under mild conditions.
Hypotheses testing, diagnostic tools and forecasting
are also proposed.
A Monte Carlo simulation is considered to evaluate the finite sample performance of the partial likelihood estimators and to study some of the proposed tests.
An empirical application is also presented and discussed.


 Keywords: double bounded time series ·  long-range dependence ·  partial likelihood · asymptotic theory ·  forecast.

MSC2000 subject classification: 62M10 · 62F12 · 62J12 · 62J99.
 
 
 
 






§ INTRODUCTION



In this work we are interested in time series whose values are restrained
to a continuous interval of the real line, say (a,b), where a<b.
Without loss of generality we consider series in the unit interval (0,1).
One typical broad case is when the time series represent rates and proportions observed over time. 
Building over the works of <cit.>, <cit.> and <cit.>,  <cit.> introduces
the class of beta autoregressive moving average models (βARMA),
which can be viewed as a specialization of the
generalized autoregressive moving average models (GARMA) <cit.>
for beta distributed variates.







Applications of the  βARMA model spam over several areas, such as medicine <cit.>, online monitoring <cit.>, neuroscience <cit.>, among many others.





Let {y_t}_t=1^∞ be a process of interest and, aiming towards prediction and the use of partial likelihood inference, let x_t-1' denote the l-dimensional vector of  (exogenous random) covariates at time t-1 and any non-random component up to time t, to be considered in the model (a possible intercept will be considered in the model separately). Let _t-1 denote the σ-field generated by the past and present (when known) explanatory variables and possibly past values of the response variable, if they are included in the model. In this framework, the σ-field _t-1 represents all the observer's knowledge about the model up to time t-1, with a possible addition of predetermined variables at time t.


Inference in  the context of βARMA process <cit.> is conducted using a conditional likelihood approach, which only allows for deterministic covariates to be introduced in the model. In this work, we adopt the more general  approach of partial likelihood, which allows for x_t-1' to contain deterministic covariates, as in the conditional likelihood approach, but also enables the inclusion of (time dependent) random covariates, as well as any type of interaction or a mixture of these. For further details on partial likelihood inference we refer the reader to <cit.>, <cit.> and <cit.>. For details on partial likelihood in time series following generalized linear models, we refer to <cit.> and references therein.


This work is concerned with an observation-driven model in which the random component follow a conditional beta distribution, parameterized as <cit.>:

    f(y_t;μ_t,ν|_t-1)=(ν)/(νμ_t)(ν(1-μ_t)) y_t^νμ_t-1(1-y_t)^ν(1-μ_t)-1,

for 0<y_t<1,
0<μ_t<1,
and
ν>0,
where (y_t|_t-1)=μ_t and (y_t|_t-1)=μ_t(1-μ_t)/1+ν.
We note that ν is a precision parameter in the sense that the greater the ν, the smaller the variance of the distribution.
The systematic component follows the usual approach of GLM with an additional dynamic term.
Let g(·) be a twice differentiable monotonic one-to-one link function for which the inverse link is of class 𝒞^2() (the class of twice continuously differentiable functions in ).
Consider the additive specification

    g(μ_t)=η_t=x_t-1'β+τ_t,

where β'=(β_1,…,β_l) are the coefficient related to the covariates
and η_t is the linear predictor.
The particular form of τ_t is discussed in several papers <cit.>.
In βARMA models <cit.>,
τ_t is assumed to follow an ARMA(p,q) process of the type

    τ_t=∑_j=1^pϕ_j(g(y_t-j)-x_t-j-1'β)+∑_k=1^qθ_kr_t-k,

where p and q and ϕ=(ϕ_1,…,ϕ_p)' and θ=(θ_1,…,θ_q)' denote the order and coefficients of the autoregressive and moving average parts of the model, respectively, and r_t denotes an error term. When x_t-1 is non-random, this is the βARMA model of <cit.>.

An advantage of specification (<ref>) is that even though the conditional mean is transformed, it is actually μ_t that is being modeled. In some applications it is common to model g(y_t) and then transform estimates  back by applying g^-1 which can be problematic (Jensen's inequality, delta method, etc.).
Observe that the time series part of the βARMA model can only accommodate short range dependence, structure that may not be enough in certain situations. In this work we propose a generalization of the βARMA model of <cit.>  by allowing τ_t to accommodate long-range dependence.






§ THE PROPOSED MODEL



The most widely applied model for time series presenting long-range dependence is the class of autoregressive fractionally integrated moving average (ARFIMA) models, introduced by
<cit.> and <cit.> <cit.>.
Recall that a process {z_t}_t=1^∞ is called an ARFIMA(p,d,q) if it is a weakly stationary solution of 

    ϕ(L)(1-L)^dz_t=θ(L)ε_t.

Here L denotes the backward shift operator L^k(z_t)=z_t-k, for k∈{1,2,…}, ε_t is an error term (usually taken as a white noise), ϕ(z) and θ(z) denote the AR and MA polynomials given respectively by

    ϕ(z)=-∑_i=0^pϕ_iz^i,    θ(z)=∑_j=0^qθ_jz^j,   ∀ z∈,

assumed, as usual, to present no common roots, where ϕ_0=-1 and θ_0=1. The fractional term (1-L)^d is defined by its binomial expansion, which, in more useful form reads

    (1-L)^-d=∑_k=0^∞π_kL^k,    π_k=Γ(k+d)/Γ(k+1)Γ(d)=∏_j=1^k j-1+d/j,  k≥1,

and π_0=1. In this work we shall assume d∈(-0.5,0.5). In this range, it can be shown that if the polynomial ϕ(z) does not have roots in the unitary disk {z∈:|z|=1}, then the ARFIMA(p,d,q) is weakly stationary.
More details on the theory of ARFIMA processes can be found in <cit.> and <cit.>.

In this work we propose a generalization of the βARMA model by allowing τ_t to follow an ARFIMA(p,d,q) process. To motivate the model, following a similar approach as <cit.>,  assume that, conditionally to _t-1, {g(y_t)-x_t-1'β}_t=1^∞ is a zero-mean stationary ARFIMA(p,d,q) process and write

    ϕ(L)(1-L)^d(g(y_t)-x_t-1'β)=θ(L)r_t 
    
       ⟹   g(y_t)=x_t-1'β+∑_j=1^pϕ_j(g(y_t-j)-x_t-j-1'β)+r_t+∑_k=1^∞ c_kr_t-k,

where we define θ_k=0, for k>q and r_t is an _t-measurable error term satisfying (r_t|_t-1)=0, for all t>1. The error term is defined in a recursive fashion in the prediction scale, that is, we consider r_t=g(y_t)-g(μ_t).

The c_k's in (<ref>)  are the coefficients obtained from Laurent's expansion of (1-z)^-dθ(z), namely

    c_0=1,       c_k=∑_i=0^min{k,q}θ_iπ_k-i,  k> 0.

Taking conditional expectation with respect to _t-1 in (<ref>), noticing that (g(y_t)-x_t-1'β|_t-1)≈τ_t,


upon substituting τ_t=g(μ_t)-x_t-1'β and adding an intercept α∈ for g(μ_t), we arrive at

    η_t=g(μ_t)=α + x_t-1'β+∑_j=1^pϕ_j(g(y_t-j)-x_t-j-1'β)+∑_k=1^∞ c_kr_t-k,

with c_k given in (<ref>). Observe that the righthand side of (<ref>) is a real number, hence η_t∈ for all t. Specification (<ref>) together with (<ref>) define the proposed βARFIMA(p,d,q) model.




§ PARAMETER ESTIMATION


In this section we shall derive the partial maximum likelihood estimator for the parameters in the proposed βARFIMA model. Let {(y_t,x_t^')}_t=1^n a sample from a βARFIMA(p,d,q) model. Let us denote the (p+q+l+3)-dimensional parameter vector by  γ'=(ν, d,α, β', ϕ',θ')  and let Ω⊆ (0,∞)× (-0.5,0.5)×^p+q+l+1 be the parameter space.
By letting

    ℓ_t(μ_t,ν)   =log(f(y_t;μ_t,ν|_t-1))
       =log((ν))-log((μ_tν))-log((ν(1-μ_t)))+
              +(μ_tν-1)log(y_t)+(ν(1-μ_t)-1)log(1-y_t),

the partial log-likelihood function is given by

    ℓ(γ)=∑_t=1^nℓ_t(μ_t,ν)

and hence, the partial maximum likelihood estimator of γ is defined as

    γ=γ∈Ωargmax{ℓ(γ)}.

In the next section we shall derive the score vector related to the maximization problem (<ref>).




 §.§ Partial score vector


To derive the partial score vector we shall need to obtain the derivative of the log-likelihood ℓ(γ) given in (<ref>) with respect to each coordinate γ_j of the parameter γ'=(ν, d,α, β', ϕ',θ').
Let ψ:(0,∞)→ be the digamma function defined as ψ(z)=d/dzlog(Γ(z)). The derivative of ℓ(γ) with respect to ν can be easily obtained as

    ∂ℓ(γ)/∂ν=∑_t=1^n [μ_t(y_t^∗-μ_t^∗)+log(1-y_t)-ψ((1-μ_t)ν)+ψ(ν)],

where

    y_t^ ∗ = log(y_t/1-y_t)     μ_t^ ∗ = ψ(μ_tν)-ψ((1-μ_t)ν).

The derivative of the log-likelihood ℓ(γ) with respect to the remaining parameters γ_j ≠ν can be computed by the general differentiation rule

    ∂ℓ(γ)/∂γ_j =
     ∑_t=1^n∂ℓ_t(μ_t,ν)/∂μ_td μ_t/d η_t∂η_t/∂γ_j.

Observe that

    ∂ℓ_t(μ_t,ν)/∂μ_t = ν[
    log(y_t/1-y_t) - ψ (μ_t ν) + ψ( ν(1-μ_t) )
    ]=ν(y_t^∗-μ_t^∗).

Since η_t=g(μ_t), it also follows that  d μ_t/d η_t = 1/g'(μ_t). Substituting these results in (<ref>), we obtain

    ∂ℓ(γ)/∂γ_j =
     ∑_t=1^nν(y_t^∗-μ_t^∗)/g^'(μ_t)∂η_t/∂γ_j .

Hence, the task of computing the derivatives of ℓ(γ) greatly simplifies to determining the derivatives ∂η_t/∂γ_j, for each coordinate γ_j ≠ν of the parameter vector γ.  To obtain the derivative with respect to d, recall that r_t=g(y_t)-g(μ_t) and notice that

    η_td =∑_k=1^∞[r_t-kc_kd+c_kr_t-kd]=∑_k=1^∞[r_t-kc_kd-c_kη_t-kd].

For π_k given in (<ref>), by using the identity Γ^'(x)=Γ(x)ψ(x), it follows that

    π_md   =1/Γ(m+1)Γ(d)^2[Γ(d)Γ(d+m)d-Γ(d+m)Γ(d)d]
       =π_m[ψ(d+m)-ψ(d)].


Hence



    η_td = ∑_k=1^∞(r_t-k∑_i=0^min{k,q}θ_iπ_k-i[ψ(d+k-i)-ψ(d)]-c_kη_t-kd).

Differentiation with respect to α yields

    ∂η_t/∂α= 1 + ∑_k=1^∞ c_k∂ r_t-k/∂α =1 - ∑_k=1^∞ c_k∂η_t-k/∂α.







Regarding parameter β_s, for s∈{1,…,l}, we have

    η_tβ_s=x_t-1(s)-∑_j=1^p ϕ_jx_t-j-1(s)-∑_k=1^∞ c_kη_t-kβ_s,

where x_m(s) denotes the (s)-th element of x_m.


The log-likelihood derivative with respect to ϕ_s, for s∈{1,…,p}, is given by

    η_tϕ_s=g(y_t-s)-x_t-s-1^'β-∑_k=1^∞ c_kη_t-kϕ_s.






For the log-likelihood derivative with respect to θ_s, with s=1,…,q, we have

    η_tθ_s=∑_k=1^∞[r_t-kc_kθ_s+c_kr_t-kθ_s].

Now, differentiating (<ref>) with respect to θ_s, we obtain c_kθ_s=π_k-sI(k≥ s) so that, from (<ref>),

    ∂ℓ(γ)/∂θ_s = ∑ _ t=1 ^ nν(y_t^∗ - μ_t^∗)/g^'(μ_t)[∑_k=s^∞π_k-sr_t-k-∑_k=1^∞ c_kη_t-kθ_s].

Finally, let M be the n× l matrix whose (t,s)-th element is given by

    M_t,s=η_tβ_s=x_t-1(s)-∑_k=1^p ϕ_kx_t-k-1(s)-∑_k=1^∞ c_kη_t-kβ_s,

P be the n× p matrix whose (t,s)-th element is given by

    P_t,s=η_tϕ_s=g(y_t-s)-x_t-s-1^'β-∑_k=1^∞ c_kη_t-kϕ_s

and Q be the n× q matrix whose (t,s)-th element is given by

    Q_t,s=η_tθ_s=∑_k=s^∞π_k-sr_t-k-∑_k=1^∞ c_kη_t-kθ_s.

Let 𝒴=(y_1^∗ - μ_1^∗,…,y_n^∗ - μ_n^∗)^', =(η_1d,…,η_nd)', a=(η_1α,…,η_nα)' and T=diag{g'(μ_1)^-1,…,g'(μ_n)^-1}.
Then, the score vector can be written in matrix form as

    U(γ)=(U_ν(γ), U_d(γ),U_α(γ),U_β(γ)^',U_ϕ(γ)^',U_θ(γ)^')^'∈^p+q+l+3,

where

    U_d(γ)   =ν'T𝒴,        
    U_α(γ)=νa'T𝒴,         
    U_β(γ)=ν M'T𝒴, 
    
    U_ϕ(γ)   =ν P'T𝒴,         
    U_θ(γ)=ν Q'T𝒴,
       U_ν(γ)   =∑_t=1^n [μ_t(y_t^∗-μ_t^∗)+log(1-y_t)-ψ((1-μ_t)ν)+ψ(ν)].


The conditional maximum likelihood is obtained by
numerically
solving the nonlinear system U(γ)=0, where 0 denotes the null vector in ^p+q+l+3. 



 §.§ Conditional information matrix


In this section we derive the Fisher conditional information matrix, which will be useful later on deriving the asymptotic properties of the partial maximum likelihood estimator for the proposed model.
For i,j∈{1,…,p+q+l+3}, it can be shown that

    ∂^2ℓ(γ)/∂γ_i ∂γ_j   = ∑_t=1^n∂/∂μ_t( ∂ℓ_t(μ_t,ν)/∂μ_td μ_t/d η_t∂η_t/∂γ_j)
    d μ_t/d η_t∂η_t/∂γ_i
       = ∑_t=1^n[ ∂^2 ℓ_t(μ_t,ν)/∂μ_t^2d μ_t/d η_t∂η_t/∂γ_j
    + ∂ℓ_t(μ_t,ν)/∂μ_t∂/∂μ_t(d μ_t/d η_t∂η_t/∂γ_j) ]
    d μ_t/d η_t∂η_t/∂γ_i .

Under the regularity conditions presented in Section <ref>, (∂ℓ_t(μ_t,ν)/∂μ_t | _t-1)=0, so that

    ( . ∂^2ℓ(γ)/∂γ_i ∂γ_j| _t-1)
       = ∑_t=1^n( . ∂^2 ℓ_t(μ_t,ν)/∂μ_t^2| _t-1)
    (d μ_t/d η_t)^2
    ∂η_t/∂γ_j∂η_t/∂γ_i.

Deriving (<ref>) with respect to μ_t twice, we obtain

    ∂^2 ℓ_t(μ_t,ν)/∂μ_t^2 = -ν^2
    {ψ' (μ_t ν) +
    ψ'
    [
    (1-μ_t)ν]
    }.

Set w_t=ν^2 {ψ' (μ_t ν) + ψ' [ (1-μ_t)ν] }. Since μ_t is _t-1-measurable, we obtain 

    ( . ∂^2ℓ(γ)/∂γ_i ∂γ_j| _t-1)
    = - ∑_t=1^nw_t/g'(μ_t)^2∂η_t/∂γ_j∂η_t/∂γ_i.


Direct differentiation of ∂ℓ(γ)/∂ν with respect to γ_i, yields

    ∂^2 ℓ(γ)/∂ν∂γ_i = ∑ _ t=1 ^ n[ (y_t^∗ - μ_t^∗) - ν∂μ_t^∗/∂ν]
    1/g'(μ_t)∂η_t/∂γ_i,

where ∂μ_t^∗/∂ν =
ψ' (μ_t ν) μ_t
-ψ'
[
(1-μ_t)ν] (1-μ_t).
Under some regularity conditions (see Section <ref>), we have (y^∗_t | _t-1)=μ_t^∗, and thus

    ( . ∂^2 ℓ(γ)/∂ν∂γ_i| _t-1)
       = - ∑_t=1^nv_t/g'(μ_t)∂η_t/∂γ_i,

where v_t = ν∂μ_t^∗/∂ν = ν{ψ' (μ_t ν) μ_t -ψ'[(1-μ_t)ν] (1-μ_t)}. The expected value of the second derivative of ℓ(γ) with respect to ν is given by

    ( . ∂^2 ℓ(γ)/∂ν^2| _t-1) =
    - ∑_t=1^nS_t,

where S_t=ψ' (μ_tν)μ_t^2+ψ'[(1-μ_t)ν](1-μ_t)^2-ψ'(ν).

Finally, let v=(v_1,…,v_n)',  W = diag{w_1/g'(μ_1)^2,…,w_n/g'(μ_n)^2},

and S=diag{S_1,…,S_n}, in matrix form we have


    [ E( . ∂^2 ℓ(γ)/∂ν^2| _t-1)=-(S),                               XXXX                                  E( . ∂^2 ℓ(γ)/∂ν∂ d | _t-1)=-v'T,;             E( . ∂^2 ℓ(γ)/∂ν∂α| _t-1)=-v'Ta,                  XXXX                                   E( . ∂^2 ℓ(γ)/∂β∂ν| _t-1)=-M'Tv,;              E( . ∂^2 ℓ(γ)/∂ϕ∂ν| _t-1)=-P'Tv,                 XXXX                                   E( . ∂^2 ℓ(γ)/∂θ∂ν| _t-1)=-Q'Tv,;                  E( . ∂^2 ℓ(γ)/∂ d^2| _t-1)=-'W,              XXXX                                  E( . ∂^2 ℓ(γ)/∂ d ∂α| _t-1)=-'Wa,;                      E( . ∂^2 ℓ(γ)/∂β∂ d  | _t-1)=-M'W,       XXXX                                 E( . ∂^2 ℓ(γ)/∂ϕ∂ d  | _t-1)=-P'W,;                       E( . ∂^2 ℓ(γ)/∂θ∂ d  | _t-1)=-Q'W,      XXXX                                    E( . ∂^2 ℓ(γ)/∂α^2| _t-1)=-(W), ]


    [          E( . ∂^2 ℓ(γ)/∂β∂α| _t-1)=-M'Wa,     XXXX                   E( . ∂^2 ℓ(γ)/∂ϕ∂α| _t-1)=-a'WP,;           E( . ∂^2 ℓ(γ)/∂θ∂α| _t-1)=-Q'Wa,    XXXX                  E( . ∂^2 ℓ(γ)/∂β∂β'| _t-1)=-M'WM,;  E( . ∂^2 ℓ(γ)/∂β∂ϕ'| _t-1)=-M'WP,            XXXX                  E( . ∂^2 ℓ(γ)/∂β∂θ'| _t-1)=-M'WQ,; E( . ∂^2 ℓ(γ)/∂ϕ∂ϕ'| _t-1)=-P'WP,             XXXX                  E( . ∂^2 ℓ(γ)/∂ϕ∂θ'| _t-1)=-P'WQ,;     E( . ∂^2 ℓ(γ)/∂θ∂θ'| _t-1)=-Q'WQ.         XXXX                                                    ]

Now, for i,j∈{1,…,p+q+l+3} let

    G_γ_i,γ_j(γ)=-E( . ∂^2 ℓ(γ)/∂γ_i ∂γ_j'| _t-1),

and G_γ_j,γ_i(γ)=G_γ_i,γ_j(γ)', for i>j.
Hence, the Fisher information matrix is the matrix G_n(γ) whose (i,j)-th element is G_γ_i,γ_j(γ).




§ ASYMPTOTIC THEORY AND HYPOTHESIS TESTING


Rigorous asymptotic theory for the maximum likelihood estimator in the context of generalized linear models for canonical link functions was first developed in <cit.> and <cit.>. For non-canonical links, the work of <cit.> was a pioneer, setting grounds for latter development of the theory. The work of <cit.> develops the theory of partial likelihood in the context of non-Gaussian and non-stationary time series.  For GARMA-like models <cit.>, a general theory for PMLE is presented in the works of <cit.>, from which we build up upon, and (at some extent) <cit.>. We remark that, since the βARMA model can be viewed as a special case of our model (when d=0 and the covariates are all non-random) the asymptotic theory presented here completes the one presented in <cit.> and <cit.>.
Let {(y_t,x_t)}_t=1^n be a sample from a βARFIMA(p,d,q) model specified by (<ref>) and (<ref>). Let  U(γ) denote the partial score vector based on the sample and γ denote a solution of U(γ)=0. Also, for j∈ let
h(t,j)=g(y_t-j)-P_{x_1',…,x_t-j-1'}(g(y_j)) where P_{x_1',…,x_t-j-1'}(g(y_j)) denotes the projection of g(y_j) into the space generated by x_1',…,x_t-j-1' and

    Z_t=(1,x_t-1',h(t,1),…,h(t,p), r_t-1,r_t-2,…)'.

To calculate the PMLE γ_n from a sample, we approximate the derivatives by truncating the infinite moving average representation (<ref>) to a point  m, initialize r_t=0 and μ_t=0 for t≤ p and calculate μ_t and r_t for t>p recursively from the data through (<ref>).


The required regularity conditions for the asymptotic existence/uniqueness, the consistency and the asymptotic normality of the partial likelihood estimator for the βARFIMA models are fundamentally the same as in <cit.>.



  §.§.§ Assumptions




  * The inverse link function g^-1 is of class 𝒞^2 and satisfies |∂ g^-1(x)/∂ x|≠0, for all x∈.

  * The parametric space Ω is an open set in ^l+q+p+3 and the true parameter γ_0 lies in Ω.

  * For each t, the covariate vector Z_t almost surely belongs to a compact set Υ⊂Ω and there exists n_0∈ such that, for all n>n_0, P(∑_t=1^nZ_tZ_t'>0)=1. Additionally, assume that g^-1(η_t) is almost surely well-defined for all Z_t∈Υ and γ∈Ω.

  * There exists a probability measure λ in Ω such that ∫_Ωzz'λ(dz) is positive definite and such that the weak convergence
    
    1/n∑_t=1^nI(Z_t-1∈ A)n→∞⟶λ(A),

holds for all λ-continuity sets A⊂Ω under (<ref>) with γ=γ_0.


Assumptions A and B guarantee that ∂^2 ℓ(γ)/∂γ'∂γ is a continuous function of γ, while Assumptions B and C imply that, for all sufficiently large n, the conditional information matrix is positive definite. Assumptions A, B, and C also assure that the model is well defined. The compactness assumption in C is mathematically convenient. It can, however, be replaced by the requirement that there exists an increasing sequence of compact sets, {Υ_n}_n=1^∞, say, such that, for sufficiently large n, x_n∈Υ_n with high probability. A probability measure satisfying the weak convergence in Assumption D also satisfies

    1/n∑_t=1^nf( Z_t-1)Pn→∞⟶∫_Υ f( z)λ(d z),

for all bounded and continuous function f:Υ→,
and implies the weak convergence of the conditional information matrix to a (non-random) positive definite matrix in the sense that there exists a matrix, which we denote by G(γ), such that

    G_n(γ)/nn→∞⟶ G(γ),  ∀γ∈Ω.

Observe that assumption D implies that G(γ_0) is positive definite and its inverse exists. Conditions C and D also imply conditions C and D in <cit.>, which, in turn, imply the asymptotic existence of a sequence of solutions for U(γ)=0. See also the discussion on <cit.> and <cit.>.


Under the assumptions A-D, the probability that a locally unique maximum partial likelihood estimator exists in a neighborhood of γ_0 tends to one. Furthermore, the estimator is consistent

    γPn→∞⟶γ_0

and asymptotically normal

    √(n)(γ-γ_0)dn→∞⟶ N_p+q+l+3(0,G(γ_0)^-1).





The proof follows the same lines as the proof of Theorem 3.1 in <cit.> by the _t-measurability of r_t and since, under the hypothesis, h(t,j)→x_t-j-1'β, almost surely. The key point is to show that the score vector is a zero mean square integrable martingale sequence with respect to an adequate filtration. Consider the filtration {_t, t∈} where _t=σ{Y_t, Y_t-1,…, x_t, x_t-1,…} and let {U_t(γ)}_t∈ denote the partial score process given by

    U_t(γ):=(∂ℓ_t(γ)/∂γ_1,…,∂ℓ_t(γ)/∂γ_l+q+p+3)^'.

First observe that the _t-1-measurability of μ_t and assumptions A and B imply that {U_t(γ)}_t∈ is integrable and adapted to the filtration {_t, t∈}. For γ_j≠ν it is straightforward to show that, under assumptions A and B, (y_t^∗|_t-1)=μ_t^∗, hence, by (<ref>), we conclude that (∂ℓ_t(μ_t,ν)/∂μ_t |_t-1)=0. When γ_j=ν, the result follows from (<ref>) since (log(1-y_t)|_t-1)=ψ((1-μ_t)ν)-ψ(ν). The rest of the proof follows the same idea as Theorem 3.1 in <cit.>, in view of Theorem 1 in <cit.> and <cit.>.
In view of the work of <cit.>, one could also, in principle, obtain similar large sample results for the partial likelihood considering non-stationary ARFIMA processes in (<ref>), under somewhat more stringent conditions.


It is widely known that, under long-range dependence, convergence rates of central limit type theorems are usually slower than the traditional √(n). Observe, however, that a time series {y_t} following a βARFIMA model is a sequence of conditionally independent, but not identically distributed random variables. In the presence of time dependent covariates, it is, in fact, non-stationary. The long-range dependence is connected to g(μ_t), which does not influence the convergence rate of the corresponding parameters in Theorem <ref>.


 Let T:^p+q+l+3→^k, for k<p+q+l+3, be a vector valued transformation such that its Jacobian J(γ) exists, is of full rank k and, as a function of γ, is continuous in an open subset of Ω. We shall consider composite hypothesis of the form

    ℋ_0:T(γ)= 0     ℋ_1:T(γ)≠ 0.

There are several ways to test the restriction (<ref>). Let γ be the unrestricted PMLE of γ

and let γ̃ be the PMLE under ℋ_0 in (<ref>).
The partial log-likelihood ratio statistic is given by

    LR=2[ℓ(γ)-ℓ(γ̃)].

The traditional Wald's statistic reads

    W=nT(γ)'[J(γ)'G^-1(γ)J(γ)]^-1T(γ),

while the Rao's score statistic is given by

    S=1/nU(γ̃)'G^-1(γ̃)U(γ̃).

Next theorem shows that the asymptotic distribution of the test statistics LR, W, and S are analogous as their counterparts under independence. The proof is completely analogous to the independent case <cit.>.

Under assumptions A-D and under the null hypothesis in (<ref>), the test statistics LR, W, and S defined above are asymptotically distributed as chi-square with k degrees of freedom.


The square root of the traditional Wald's statistic (often called z statistics)
is particularly convenient to test individual parameters <cit.>.
Considering the hypothesis
ℋ_0:γ_j= γ_0 vs. ℋ_1:γ_j ≠γ_0,
the z statistic is given by

    z= γ_j - γ_0/se(γ_j),

where se(γ_j) is the square root of the j-th diagonal element of G_n(γ)^-1.
Under ℋ_0, the limiting distribution of z is standard normal.

As an example the transformation T(γ)=d coupled with any of the above test statistics can be used to test

    ℋ_0: d=0     ℋ_1: d≠0,

which is equivalent to test the presence of long-range dependence in the systematic component of the model. In other words, it can be applied to decide whether a βARFIMA or a βARMA is suitable to the data, with rejection of the null hypothesis favouring the βARFIMA model.



§ DIAGNOSTIC AND PREDICTION


Diagnostics in the context of  βARFIMA models follow the usual procedures of GLM theory with some adaptations. For a general goodness of fit testing we consider the so-called deviance statistic.
The deviance D is defined as twice the difference between the conditional log-likelihood of the saturated model (for which μ̃_t=y_t, i.e., a model with as many parameters as observations) and the fitted model, that is

    D=2( ℓ̃ - ℓ ),    ℓ  =  ∑_t=1^nℓ_t(μ_t,ν)  ℓ̃  =  ∑_t=1^nℓ_t(y_t,ν).


If the fitted model is correct,
the test statistic D is approximately distributed as chi-squared with n-(p+q+l+3) degrees of freedom <cit.>.






Model selection among several competing models may be based on the usual information criteria. The Akaike information criterion (AIC) <cit.> is given by

    AIC = -2 ℓ + 2 (p+q+l+3).

Usual information criteria aims at estimating the expected partial log-likelihood and applying a penalty proportional to the number of parameters in the model (in the AIC case, 2(p+q+l+3)) for the maximized partial log-likelihood function. If the penalty term 2(p+q+l+3) in (<ref>) is replaced by log(n)(p+q+l+3), we obtain the Schwarz information criterion (BIC) <cit.>; if it
is replaced by log(log(n))(p+q+l+3) instead, the <cit.> criterion (HQ) is obtained.



Residual analysis is an important step in verifying whether the estimated model provides a good fit to the data.
Since the proposed model is an extension of the beta regression model <cit.>,
the residual analysis applied to the former can be also applied to the βARFIMA models <cit.>.
At the outset, we define the following standardized residual:


    r_t^1 = y_t-μ_t/√( Var(y_t)) =
    y_t - μ_t/√(μ_t(1-μ_t)/(1+ν)),    t=1,…, n.

A more sophisticated residual is the standardized weighted residual, introduced by <cit.>, which is given by

    r_t^w = y_t^∗-μ_t^∗/√( Var(y_t^∗)) =
    y_t^∗-μ_t^∗/√(ψ'(μ_t ν) + ψ'[(1-μ_t) ν] ),    t=1,…, n.

The authors have shown that Var(y_t^∗)= ψ'(μ_t ν) + ψ'[(1-μ_t) ν]. See also <cit.>.


When the fitted model is correct, these residuals are well approximated by the standard normal distribution.

When the model is correctly specified, the residuals should display white noise behavior, i.e., they should follow a zero mean and constant variance uncorrelated process <cit.>. A visual inspection of the residual plot is an indispensable tool for a first step residual check <cit.>.
Let r_1^(·), …, r_n^(·) be any type of residual obtained from the fitted model.
The usual estimate for the residual autocorrelation function (ACF) is

    ρ(h)= ∑_t=1^n-h(r_t^(·)-r^(·))(r_t+h^(·)-r^(·))/∑_t=1^n-h(r_t^(·)-r^(·))^2,    h=0,1,…,

where r^(·) = 1/n∑_t=1^nr_t^(·). When i>1 and n is sufficiently large, the distribution of ρ(i) is approximately normal with zero mean and variance 1/n <cit.>.
Hence, one can apply the usual ± 1.96/√(n) as (95%) confidence bands in ACF plots as a first visual inspection tool for white noise behavior <cit.>. Since these bounds are usually conservative in finite samples (tighter than they should be), a Ljung-Box test can also be applied. In that case, for large enough n, the test statistics will follow its usual distribution, but see the discussion in <cit.>.


Applying the partial maximum likelihood estimator  in (<ref>), we can obtain h steps ahead predicted values for the observed response y_t, which we denote by  y_n+h=y_n(h),


    y_n(h)= g^-1(α + x_t-1'β
    + ∑_j=1^p ϕ_j[g(y_n+h-j)-x_n+h-j-1'β]
    + ∑_k=1^m c_k [r_n+h-k]
    ),


where m is the (user chosen) truncation point for the MA(∞) representation in (<ref>), c_k is the quantity in (<ref>) evaluated at the PMLE estimates θ and d,


    [g(y_n+h-j)-x_n+h-j-1'β]
        =
    
    g(y_n(h-j))-x_n+h-j-1'β,   if j< h,
    
    g(y_n+h-j)-x_n+h-j-1'β,    if j≥ h,


and


    [r_t]    =
    
    0,   if t≤ 0 or t≥ n+h,
    
    g(y_t) - g(μ_t),    if 1≤ t≤ n,
    
    g(y_n(t-n)) - g(μ_t),    if n+1≤ t ≤ n+h-1.


Finally, in the absence of covariates or when each of the covariates forms a stationary sequence with absolutely summable autocorrelation function, identification of long-range dependence in the conditional mean can be done by using the covariance decay of either y_t, g(y_t) or cumulative average as a proxy. A slow ACF decay in any of these sequences indicates the presence of long-range dependence in the conditional mean. In the presence of time varying covariates with either, non-stationary, non absolutely summable ACF  or deterministic behavior, the diagnostic can only be done after dealing with the non-stationarity in the series.




§ MONTE CARLO SIMULATION



In this section we  present a Monte Carlo simulation study to assess the finite sample properties of the PMLE for βARFIMA models
as well as the LR and Wald's z tests for the presence of long-range dependence. Observe that the z statistics is obtained from the information matrix under the alternative hypothesis, while the LR statistics is directly obtained from the log-likelihood, which favor their use in detriment of the Rao's Score test, which requires matrix inversion and evaluation of the score vector under both hypothesis.
We simulate 1,000 replicates of a βARFIMA(1,d,1) model restricted to the interval (0,1),
with
ϕ_1=0.2, θ_1=-0.3, two values of ν∈{40,120}, d∈{0.15,0.30,0.45} and sample sizes n∈{1000,3000,5000}. We apply the logit as link function and no covariates were included in the simulations.

Given the vector of parameters γ, to generate a size n sample from the specified βARFIMA(p,d,q) process restricted to the interval (0,1), let m>p denote the cutoff point for the infinite sum in (<ref>). We start the algorithm by setting r_t=0 for all t≤ m and μ_t=g^-1(α), for t=1,…,m. Second step: for t=m+1, we obtain η_t through (<ref>), then we set μ_t=g^-1(η_t) and update r_t=g(y_t)-g(μ_t). Finally, y_t is generated from (<ref>), using any adequate method (such as the inversion method). We iterate the second step for t=m+1,…,n_0+n, where n_0>m denotes the size of a possible burn in. The desired sample is y_n_0+1,…,y_n_0+n. If needed, the sample can be rescaled to (a,b) through ỹ_t=a+(b-a)y_t. We have also performed a pilot simulation study (not shown) to determine the influence of the cutoff point m in parameter estimation. We found that for m≥50, it has negligible impact on the estimated values and that a good compromise between computational speed and accuracy is m=100 (used here). For practical purposes, where only a handful of series are analyzed, m=200 seems a good choice.  All routines were implemented by the authors and are available in R language <cit.> upon request. The code for the main tasks of computing the partial score vector and the information matrix were written in FORTRAN 90 by the authors and called from within R (βARMA models can also be fitted). Optimization is performed by using the so-called L-BFGS-S algorithm <cit.>, which was also implemented in FORTRAN 90 language based on <cit.> and applied without any parameter constraint. We use analytical derivatives in the optimization procedure, given in Section <ref>. The iterative optimization algorithm requires initialization. The starting values of the constant (α),
the covariate parameter (β)
and the autoregressive (ϕ) parameters
were obtained from a linear regression with response
 Y=(g(y_m+1), g(y_m+2), …, g(y_n))^' on the design matrix X

    X=
    [          1       x_m1       x_m2          ⋯       x_mr     g(y_m)   g(y_m-1)          ⋯ g(y_m-p+1);          1   x_(m+1)1   x_(m+1)2          ⋯   x_(m+1)r   g(y_m+1)     g(y_m)          ⋯ g(y_m-p+2);          ⋮          ⋮          ⋮          ⋱          ⋮          ⋮          ⋮          ⋱          ⋮;          1       x_n1       x_n2          ⋯       x_nr   g(y_n-1)   g(y_n-2)          ⋯   g(y_n-p);            ],

For the parameter θ, the starting values are set to zero and d is started as 0.001. We apply analytic derivatives where the ones obtained through iteration are initialized with zero for non-observed values. 









Tables <ref> and <ref> present the simulation results for point estimates. Performance statistics presented are the mean, percentage relative bias (RB%), variance (Var) and mean square error (MSE).
The percentage relative bias is defined as the ratio between the bias and the true parameter value times 100.
For n=1,000, there is a small bias for parameter d, which is expected since, in the context of long-range dependent processes, it is quite common the presence of bias for smaller sample sizes <cit.>.
Overall, the results in Table <ref>, for all d's show somewhat smaller bias for the parameters α,ν,d and considerably higher bias for the other estimates. As expected, as n increases, the bias in the estimates decrease (except for ν, but the difference is so small that it can be considered negligible) and so are the variance and MSE, which is a reflection of the PLME's consistency.







We also evaluate the performance of the z and LR statistics for testing the null hypothesis ℋ_0:   d=0
against two-sided alternative hypothesis.
For this purpose we consider three nominal levels: 1%, 5% and 10%, and the same scenarios described above.
 

 Table <ref> presents the null rejection rates of the two different tests. The figures in this table
clearly show that the test based on z statistics is considerably oversized (liberal) in smaller samples.
In the other hand, the LR test presents the best performer,  being much less distorted than z test.
The LR test's null rejection rates are closer to the nominal levels than Wald's.



We also present the non-null rejection rates, i.e., their estimated power. The results are presented in Table <ref>.
As expected, the tests become more powerful as d moves away from zero
and as n increases.
We also notice that the z test is more powerful than the LR test. However, the z test is considerably oversized and this can be an unfair comparison.
Therefore, we conclude that the LR test is more reliable to test the presence of long-range dependence in βARFIMA model than the z test.






§ REAL DATA APPLICATION


The relative air humidity (or simply relative humidity, abbreviated RH) is an important meteorological characteristic to public health, irrigation scheduling design, and hydrological studies. Low RH is known to cause health problems, such as allergies, asthma attacks, dehydration, nasal bleeding, among others, while high RH besides causing respiratory problems, is responsible for the increase in precipitation which, in excess, can cause serious consequences, such as flooding in urban areas, landslides, damages to agriculture, etc.


To exemplify the usefulness of the proposed βARFIMA model, we present an application to the daily average of the RH in Manaus, the Amazônia State capital in Brazil, from 01/01/2010 to 05/29/2017, yielding a sample size of n=2,704. The data from 05/30/2016 to 05/29/2017 is reserved to measure the out-of-sample forecast performance of the presented models. The data can be freely obtained from the Instituto Nacional de Meteorologia's (INMET - Brazilian National Institute of Meteorological Research) website (<http://www.inmet.gov.br>). The particular station from where the data was collected is situated at longitude 3^∘06' south, latitude 60^∘ west in Manaus.

Figure <ref> presents the RH time series and its autocorrelation function and partial autocorrelation function (PACF). The time series plot reveals a very distinctive seasonality, which we shall incorporate into the  βARFIMA model as a covariate. We define this covariate as the seasonal part of a Holt-Winters' decomposition (additive) of the time series <cit.>. This decomposition is also useful for out-of-sample forecasting as future values for the covariates can be trivially obtained from it. To fit the model, we use m=200, the logit as link function  and the diagnostics are based on the standardized residual defined in  (<ref>), while p-values are obtained from the Wald's z test (<ref>).  To select a model to the data, we systematically try different order βARFIMA(p,d,q) models and select the one whose parameter are all significant and whose residual does not reject the null hypothesis (using 20 lags) in the Ljung-Box test. All tests are conducted at 5% significance level. For comparison purposes, we also fit an additive Holt-Winters and a βARMA model. To fit and select the βARMA model we follow a similar approach as the βARFIMA. The routines in R and data used in this section are available upon request.


Based on the criteria explained above, we have selected a βARFIMA(0,d,1) model for the relative humidity data. Table <ref> presents the fitted βARFIMA model along with some diagnostics.





Proceeding similarly as in the βARFIMA case, we have selected a βARMA(1,2) for the data set. Table <ref> presents the fitted model along with some diagnostics. Simpler models did present all significative coefficients, but failed the Ljung-Box test. 




It is well known in the literature that a long-range dependent process can be well approximated by an ARMA process for which the roots of the autoregressive polynomial are close to the unit circle <cit.>. It is very interesting to notice that this is also reflected in the present case. Observe that the values of β_1 and ν on the fitted βARMA and βARFIMA  are close but the AR part of the fitted βARMA model present a root very close to the unitary circle (≈1.03). In terms of model selection,  all goodness of fit criteria (AIC, BIC and log-likelihood) suggest the βARMA as the best model for the data, but the difference is almost imperceptible. Also notice that, both Wald's z test and the LR test point to a significant long-range dependence parameter d.  

Table <ref> present the fitted additive Holt-Winters model. As expected, the trend coefficient is zero as the data presents no trend. For this model a Ljung-Box test shows that there is still serial dependence in the residuals (defined as observed minus fitted values) from the Holt-Winters model (p-value < 2.2×10^-16). The Holt-Winters is a predictive model so there are no diagnostics for it.





 We also present an in-sample and out-of-sample  forecasting study based on the fitted models. As mentioned before, we have reserved the data from 05/30/2016 to 05/29/2017 (365 observations) to compare with the out-of-sample forecasts obtained from the models.

The in-sample, 365 steps ahead out-of-sample forecasts for the fitted models, as well as the reserved data  are presented in Figure <ref> while Table <ref> presents some forecasting diagnostics, namely, the root mean squared error (RMSE), the mean absolute error (MAE), and  mean absolute percentage error (MAPE). The out-of-sample diagnostics were obtained from the 365 reserved values, compared to 365 step-ahead forecasts. 

In terms of in-sample forecast, the plots indicates that both, the βARFIMA and βARMA models, successfully captured the seasonal component in the data. Also the results presented in Table <ref> show that the βARMA present slightly better in-sample forecast diagnostics compared to the βARFIMA model, which is not surprising given that the fitted βARMA model presents more parameters. The in-sample forecast for Holt-Winters model seems visually poorer and this is reflected in the diagnostics as well.

 Out-of-sample results, however, present a totally different picture. For the βARFIMA the out-of-sample forecast seems to predict well the data behavior and the diagnostics are quite good. The out-of-sample forecasts for the βARMA are clearly off. This behavior is expected because the AR polynomial in the fitted model present a near-unit root, which induces a near-integrated process behavior in the model's conditional mean.  Hence, even though the βARMA presents a slightly better in-sample forecast for the data, slightly better goodness-of-fit measures, the model fails in producing adequate forecasts for the data due to the evidence of long-range dependence in the processes' conditional mean, which is balanced by a near unit root in the AR polynomial. Finally, the Holt-Winters model is capable of producing meaningful out-of-sample forecast, but they are overall poorer when compared to the βARFIMA's.











§ CONCLUSION AND FINAL REMARKS


In this work we introduce and study a dynamic time series regression model for bounded continuous random variables observed over time. The proposed model accommodates regressors through a GLM-type structure  and a long-range dependent time series structure. The proposed model generalizes the βARMA model of <cit.> by allowing the time series part of the model to present long-range dependence. The model also allows for covariates which can be random, non-random (predetermined) and even time dependent in nature. This flexibility is due to the adopted partial maximum likelihood approach for parameter estimation. Besides introducing the concept in the model's framework, we also explicitly derive the associated score vector and conditional information matrix.

The paper also presents the asymptotic theory for the proposed partial maximum likelihood estimator. We show that the estimator exists and is asymptotically unique, consistent and normally distributed. Based on the asymptotic theory, we provide hypothesis testing, diagnostics and forecast tools for the proposed model.

A Monte Carlo simulation study is presented and shows the proposed partial maximum likelihood estimator finite sample performance as well as the likelihood ratio and the Wald's z test's. The simulation show an overall good point estimation performance of the PMLE. Regarding the tests, we found that the Wald's z test is considerably oversized, while the LR test performs well in all simulated scenarios. 

Finally, an application to data regarding daily relative humidity in Manaus, Brazil is presented. We compare the βARFIMA, βARMA and Holt-Winters models in terms of goodness-of-fit measures as well as in-sample and out-of-sample forecasts. The results show that the βARFIMA model was capable of capturing the data dynamics, including seasonality.  Regarding  out-of-sample forecast, the βARFIMA was again the best one in terms of commonly used accuracy measures. The fitted βARMA  presented a near-unit root in the autoregressive polynomial with severe out-of-sample forecast implications. The application presents a scenario where the βARFIMA model is adequate and yields useful forecasts while the βARMA, although suitable for the data, fails to produce meaningful forecast due to the evidence of long-range dependence in the data's conditional mean. 




§ ACKNOWLEDGEMENTS


We gratefully acknowledge partial financial support from CNPq and FAPERGS, Brazil.
We are also grateful to two anonymous referees whose comments and suggestions contributed to improve the paper's quality.

elsarticle-harv




