carfima {carfima}

Fitting a CARFIMA(p, H, q) model via frequentist or Bayesian machinery

Description

A general-order CARFIMA(p, H, q) model for p>q is

Y_t^{(p)} -\alpha_p Y_t^{(p-1)} -\cdots- \alpha_1 Y_t = \sigma(B_{t, H}^{(1)}+\beta_1B_{t, H}^{(2)}+\cdots+\beta_q B_{t, H}^{(q+1)}),

where B_{t, H} = B_t^H is the standard fractional Brownian motion, H is the Hurst parameter, and the superscript (j) indicates j-fold differentiation with respect to t; see Equation (1) of Tsai and Chan (2005) for details. The model has p+q+2 unknown model parameters; p \alpha_j's, q \beta_j's, H, and \sigma. Also, the model can account for heteroscedastic measurement errors, if the information about measurement error standard deviations is known.

The function carfima fits the model, producing either their maximum likelihood estimates (MLEs) with their asymptotic uncertainties or their posterior samples according to its argument, method. The MLEs except \sigma are obtained from a profile likelihood by a global optimizer called the differential evolution algorithm on restricted ranges, i.e., (-0.99, -0.01) for each \alpha_j, (0, 100) for each \beta_j, and (0.51, 0.99) for H; the MLE of \sigma is then deterministically computed. The corresponding asymptotic uncertainties are based on a numerical hessian matrix calculation at the MLEs (see function hessian in pracma). It also computes the Akaike Information Criterion (AIC) that is -2(log likelihood -p-q-2). The function carfima becomes numerically unstable if p>2, and thus it may produce numerical errors.

The Bayesian approach uses independent prior distributions for the unknown model parameters; a Uniform(-0.9999, -0.0001) prior for each \alpha_j, a Uniform(0, 100) prior for each \beta_j, a Uniform(0.5001, 0.9999) prior for H for long memory process, and finally an inverse-Gamma(shape = 2.01, scale = 10^3) prior for \sigma^2. Posterior propriety holds because the prior distributions are jointly proper. It also adopts appropriate proposal density functions; a truncated Normal(current state, proposal scale) between -0.9999 and -0.0001 for each \alpha_j, a truncated Normal(current state, proposal scale) between 0 and 100 for each \beta_j, a truncated Normal(current state, proposal scale) between 0.5001 and 0.9999 for H, and fianlly a Normal(log(current state), proposal scale on a log scale) for \sigma^2, i.e., \sigma^2 is updated on a log scale. We sample the full posterior using Metropolis-Hastings within Gibbs sampler. It also adopts adaptive Markov chain Monte Carlo (MCMC) that updates the proposal scales every 100 iterations; if the acceptance rate of the most recent 100 proposals (at the end of the ith 100 iterations) smaller than 0.3, then it multiplies \exp(-\min(0.1, 1/\sqrt{i})) by the current proposal scale; if it is larger than 0.3, then it multiplies \exp(\min(0.1, 1/\sqrt{i})) by the current proposal scale. The resulting Markov chain with this adaptive scheme converge to the stationary distribution because the adjustment factor, \exp(\pm\min(0.1, 1/\sqrt{i})), approaches unity as i goes to infinity, satisfying the diminishing adaptation condition. The function carfima becomes numerically unstable if p>2, and thus it may produce numerical errors. The output returns the AIC for which we evaluate the log likelihood at the posterior medians of the unknown model parameters.

Usage

carfima(Y, time, measure.error, ar.p, ma.q, method = "mle", 
               bayes.param.ini, bayes.param.scale, bayes.n.warm, bayes.n.sample)

Arguments

Y

A vector of length k for the observed data.

time

A vector of length k for the observation times.

measure.error

(Optional) A vector for the k measurement error standard deviations, if such information is available (especially for astronomical applications). A vector of zeros is automatically assigned, if nothing is specified.

ar.p

A positive integer for the order of the AR model. ar.p must be greater than ma.q. If ar.p is greater than 2, numerical errors may occur.

ma.q

A non-negative integer for the order of the MA model. ma.q must be smaller than ar.p.

method

Either "mle" or "bayes". Method "mle" conducts the MLE-based inference, producing MLEs and asymptotic uncertainties of the model parameters. Method "bayes" draws posterior samples of the model parameters.

bayes.param.ini

Only if method is "bayes". A vector of length p+q+2 for the initial values of p \alpha_j's, q \beta_j's, H, and \sigma to implement a Markov chain Monte Carlo method (Metropolis-Hastings within Gibbs sampler). When a CARFIMA(2, H, 1) model is fitted, for example, users should set five initial values of \alpha_1, \alpha_2, \beta_1, H, and \sigma.

bayes.param.scale

Only if method is "bayes". A vector of length p+q+2 for jumping (proposal) scales of the Metropolis-Hastings steps. Each number determines how far a Metropolis-Hastings step reaches out in each parameter space. Note that the last number of this vector is the jumping scale to update \sigma^2 on a log scale. The adaptive MCMC automatically adjusts these jumping scales during the run.

bayes.n.warm

Only if method is "bayes". A scalar for the number of burn-ins, i.e., the number of the first few iterations to be discarded to remove the effect of initial values.

bayes.n.sample

Only if method is "bayes". A scalar for the number of posterior samples for each parameter.

Details

The function carfiam produces MLEs, their asymptotic uncertainties, and AIC if method is "mle". It produces the posterior samples of the model parameters, acceptance rates, and AIC if method is "bayes".

Value

The outcome of carfima is composed of:

mle

If method is "mle". Maximum likelihood estimates of the model parameters, p \alpha_j's, q \beta_j's, H, and \sigma.

se

If method is "mle". Asymptotic uncertainties (standard errors) of the MLEs.

param

If method is "bayes". An m by (p+q+2) matrix where m is the number of posterior draws (bayes.n.sample) and the columns correspond to parameters, p \alpha_j's, q \beta_j's, H, and \sigma.

accept

If method is "bayes". A vector of length p+q+2 for the acceptance rates of the Metropolis-Hastings steps.

AIC

For both methods. Akaike Information Criterion, -2(log likelihood -p-q-2). The log likelihood is evaluated at the MLEs if method is "mle" and at the posterior medians of the unknown model parameters if method is "bayes".

fitted.values

For both methods. A vector of length k for the values of E(Y_{t_i}\vert Y_{<t_i}), i=1, 2, \ldots, k, where k is the number of observations and Y_{<t_i} represents all data observed before t_i. Note that E(Y_{t_1}\vert Y_{<t_1})=0. MLEs of the model parameters are used to compute E(Y_{t_i}\vert Y_{<t_i})'s if method is "mle" and posterior medians of the model parameters are used if method is "bayes".

Author(s)

Hyungsuk Tak, Henghsiu Tsai, Kisung You

References

H. Tsai and K.S. Chan (2005) "Maximum Likelihood Estimation of Linear Continuous Time Long Memory Processes with Discrete Time Data," Journal of the Royal Statistical Society (Series B), 67 (5), 703-716. DOI: 10.1111/j.1467-9868.2005.00522.x

H. Tsai and K.S. Chan (2000) "A Note on the Covariance Structure of a Continuous-time ARMA Process," Statistica Sinica, 10, 989-998.
Link: http://www3.stat.sinica.edu.tw/statistica/j10n3/j10n317/j10n317.htm

Examples

  ##### Irregularly spaced observation time generation.

  length.time <- 30
  time.temp <- rexp(length.time, rate = 2)
  time <- rep(NA, length.time + 1)
  time[1] <- 0
  for (i in 2 : (length.time + 1)) {
    time[i] <- time[i - 1] + time.temp[i - 1]
  }
  time <- time[-1]

  ##### Data genration for CARFIMA(1, H, 0) based on the observation times. 

  parameter <- c(-0.4, 0.8, 0.2) 
  # AR parameter alpha = -0.4
  # Hurst parameter = 0.8
  # Process uncertainty (standard deviation) sigma = 0.2

  me.sd <- rep(0.05, length.time)
  # Known measurement error standard deviations 0.05 for all observations
  # If not known, remove the argument "measure.error = me.sd" in the following codes,
  # so that the default values (zero) are automatically assigned.

  y <- carfima.sim(parameter = parameter, time = time, 
                   measure.error = me.sd, ar.p = 1, ma.q = 0)  

  ##### Fitting the CARFIMA(1, H, 0) model on the simulated data for MLEs.
  
  res <- carfima(Y = y, time = time, measure.error = me.sd, 
                 method = "mle", ar.p = 1, ma.q = 0)
  
  # It takes a long time due to the differential evolution algorithm (global optimizer).
  # res$mle; res$se; res$AIC; res$fitted.values

  ##### Fitting the CARFIMA(1, H, 0) model on the simulated data for Bayesian inference.
  
  res <- carfima(Y = y, time = time, measure.error = me.sd, 
                 method = "bayes", ar.p = 1, ma.q = 0, 
                 bayes.param.ini = parameter, 
                 bayes.param.scale = c(rep(0.2, length(parameter))), 
                 bayes.n.warm = 100, bayes.n.sample = 1000)
  
  # It takes a long time because the likelihood evaluation is computationally heavy.
  # The last number of bayes.param.scale is to update sigma2 (not sigma) on a log scale.
  # hist(res$param[, 1]); res$accept; res$AIC; res$fitted.values

  ##### Computing the log likelihood of the CARFIMA(1, H, 0) model given the parameters.
  loglik <- carfima.loglik(Y = y, time = time, ar.p = 1, ma.q = 0,
                           measure.error = me.sd,
                           parameter = parameter, fitted = FALSE)