Whittle estimator to Locally Stationary Time Series

Description

This function computes Whittle estimator to LS-ARMA and LS-ARFIMA models.

Usage

LS.whittle(series, start, order = c(p = 0, q = 0), ar.order = NULL,
           ma.order = NULL, sd.order = NULL, d.order = NULL, include.d = FALSE, 
           N = NULL, S = NULL, include.taper = TRUE, control = list(),
           lower = -Inf, upper = Inf, m = NULL, n.ahead = 0)

Arguments

Details

This function estimates the parameters in models: LS-ARMA

Φ(t/T, , B), Y_{t, T} = Θ(t/T,, B),σ(t/T), varepsilon_t

and LS-ARFIMA

Φ(t/T, , B), Y_{t, T} = Θ(t/T,, B), (1-B)^{-d(t/T)}, σ(t/T), varepsilon_t,

with infinite moving average expansion

Y_{t, T} = σ(t/T), ∑_{j=0}^{∞} ψ(t/T),varepsilon_t,

for t = 1,…, T, where for u = t/T in [0,1], Φ(u,B)=1+φ_1(u)B +cdots+φ_p(u)B^p is an autoregressive polynomial, Θ(u, B) = 1 + θ_1(u)B + cdots + θ_q(u)B^q is a moving average polynomial, d(u) is a long-memory parameter, σ(u) is a noise scale factor and {varepsilon_t } is a Gaussian white noise sequence with zero mean and unit variance. This class of models extends the well-known ARMA and ARFIMA process, which is obtained when the components Φ(u, B), Θ(u, B), d(u) and σ(u) do not depend on u.

The evolution of these models can be specified in terms of a general class of functions. For example, let {g_j(u)}, j = 1, 2, …, be a basis for a space of smoothly varying functions and let d_{θ}(u) be the time-varying long-memory parameter in model LS-ARFIMA. Then we could write d_{θ}(u) in terms of the basis {g_j(u) = u^j} as follows d_{θ}(u) = ∑_{j=0}^{k} α_j,g_j(u) for unknown values of k and θ = (α_0,,α_1,,…, ,α_k)^{prime}. In this situation, estimating θ involves determining k and estimating the coefficients α_0,,α_1,,…, ,α_k.

LS.whittle optimizes LS.whittle.loglik as objective function using nlminb function, for both LS-ARMA (include.d=FALSE) and LS-ARFIMA (include.d=TRUE) models. Also computes Kalman filter with LS.kalman and this values are given in var.coef in the output.

Value

A list with the following components:

Author(s)

Ricardo Olea <raolea@uc.cl>

See Also

nlminb for optimization. LS.kalman for Kalman filter.

Examples

## Require "rdatamarket"
library(rdatamarket)
malleco = dmlist("22tn")
y = malleco$Value

# Analysis by blocks of phi and sigma parameters
T. = length(y)
N = 200
S = 100
M = trunc((T. - N)/S + 1)
table = c()
for(j in 1:M){
x = y[(1 + S*(j-1)):(N + S*(j-1))]
table = rbind(table,nlminb(start = c(0.65, 0.15), N = N, objective = LS.whittle.loglik,
                           series = x, order = c(p = 1, q = 0))$par)
}
u = (N/2 + S*(1:M-1))/T.
table = as.data.frame(cbind(u, table))
colnames(table) = c("u", "phi", "sigma")

spar = 0.6
## start parameters
phi = smooth.spline(table$phi, spar = spar)$y
fit.1 = nls(phi ~ a0+a1*u, start = list(a0 = 0.65, a1 = 0.00))
sigma = smooth.spline(table$sigma, spar = spar)$y
fit.2 = nls(sigma ~ b0+b1*u, start = list(b0 = 0.65, b1 = 0.00))

fit = LS.whittle(series = y, start = c(coef(fit.1), coef(fit.2)), order = c(p = 1, q = 0),
                 ar.order = c(1), sd.order = 1, N = 180, n.ahead = 10)

names(fit)
fit$coef
fit$loglik

Results

R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(LSTS)
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/LSTS/LS.whittle.Rd_%03d_medium.png", width=480, height=480)
> ### Name: LS.whittle
> ### Title: Whittle estimator to Locally Stationary Time Series
> ### Aliases: LS.whittle whittle estimator
> ### Keywords: whittle loglik estimator timeseries
> 
> ### ** Examples
> 
> ## Require "rdatamarket"
> library(rdatamarket)
Loading required package: zoo

Attaching package: 'zoo'

The following objects are masked from 'package:base':

    as.Date, as.Date.numeric

> malleco = dmlist("22tn")
> y = malleco$Value
> 
> # Analysis by blocks of phi and sigma parameters
> T. = length(y)
> N = 200
> S = 100
> M = trunc((T. - N)/S + 1)
> table = c()
> for(j in 1:M){
+ x = y[(1 + S*(j-1)):(N + S*(j-1))]
+ table = rbind(table,nlminb(start = c(0.65, 0.15), N = N, objective = LS.whittle.loglik,
+                            series = x, order = c(p = 1, q = 0))$par)
+ }
> u = (N/2 + S*(1:M-1))/T.
> table = as.data.frame(cbind(u, table))
> colnames(table) = c("u", "phi", "sigma")
> 
> spar = 0.6
> ## start parameters
> phi = smooth.spline(table$phi, spar = spar)$y
> fit.1 = nls(phi ~ a0+a1*u, start = list(a0 = 0.65, a1 = 0.00))
> sigma = smooth.spline(table$sigma, spar = spar)$y
> fit.2 = nls(sigma ~ b0+b1*u, start = list(b0 = 0.65, b1 = 0.00))
> 
> fit = LS.whittle(series = y, start = c(coef(fit.1), coef(fit.2)), order = c(p = 1, q = 0),
+                  ar.order = c(1), sd.order = 1, N = 180, n.ahead = 10)
> 
> names(fit)
 [1] "coef"          "var.coef"      "loglik"        "aic"          
 [5] "series"        "residuals"     "fitted.values" "pred"         
 [9] "se"            "model"        
> fit$coef
         a0          a1          b0          b1 
 0.50055762  0.25769176  0.11295694 -0.01224227 
> fit$loglik
[1] 2.662094
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>