AR1: Autoregressive Process with Order-1 Family Function

Description Usage Arguments Details Value Warning Note Author(s) References See Also Examples

View source: R/family.ts.R

Description

Maximum likelihood estimation of the three-parameter AR-1 model

Usage

1
2
3
4
5
AR1(ldrift = "identitylink", lsd  = "loglink", lvar = "loglink", lrho = "rhobitlink",
    idrift  = NULL, isd  = NULL, ivar = NULL, irho = NULL, imethod = 1,
    ishrinkage = 0.95, type.likelihood = c("exact", "conditional"),
    type.EIM  = c("exact", "approximate"), var.arg = FALSE, nodrift = FALSE,
    print.EIM = FALSE, zero = c(if (var.arg) "var" else "sd", "rho"))

Arguments

ldrift, lsd, lvar, lrho

Link functions applied to the scaled mean, standard deviation or variance, and correlation parameters. The parameter drift is known as the drift, and it is a scaled mean. See Links for more choices.

idrift, isd, ivar, irho

Optional initial values for the parameters. If failure to converge occurs then try different values and monitor convergence by using trace = TRUE. For a S-column response, these arguments can be of length S, and they are recycled by the columns first. A value NULL means an initial value for each response is computed internally.

ishrinkage, imethod, zero

See CommonVGAMffArguments for more information. The default for zero assumes there is a drift parameter to be estimated (the default for that argument), so if a drift parameter is suppressed and there are covariates, then zero will need to be assigned the value 1 or 2 or NULL.

var.arg

Same meaning as uninormal.

nodrift

Logical, for determining whether to estimate the drift parameter. The default is to estimate it. If TRUE, the drift parameter is set to 0 and not estimated.

type.EIM

What type of expected information matrix (EIM) is used in Fisher scoring. By default, this family function calls AR1EIM, which recursively computes the exact EIM for the AR process with Gaussian white noise. See Porat and Friedlander (1986) for further details on the exact EIM.

If type.EIM = "approximate" then approximate expression for the EIM of Autoregressive processes is used; this approach holds when the number of observations is large enough. Succinct details about the approximate EIM are delineated at Porat and Friedlander (1987).

print.EIM

Logical. If TRUE, then the first few EIMs are printed. Here, the result shown is the sum of each EIM.

type.likelihood

What type of likelihood function is maximized. The first choice (default) is the sum of the marginal likelihood and the conditional likelihood. Choosing the conditional likelihood means that the first observation is effectively ignored (this is handled internally by setting the value of the first prior weight to be some small positive number, e.g., 1.0e-6). See the note below.

Details

The AR-1 model implemented here has

Y(1) ~ N(mu, sigma^2 / (1-rho^2),

and

Y(i) = mu^* + rho * Y(i-1) + e(i)

where the e(i) are i.i.d. Normal(0, sd = sigma) random variates.

Here are a few notes: (1). A test for weak stationarity might be to verify whether 1/rho lies outside the unit circle. (2). The mean of all the Y(i) is mu^* / (1-rho) and these are returned as the fitted values. (3). The correlation of all the Y(i) with Y(i-1) is rho. (4). The default link function ensures that -1 < rho < 1.

Value

An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm, and vgam.

Warning

Monitoring convergence is urged, i.e., set trace = TRUE.

Moreover, if the exact EIMs are used, set print.EIM = TRUE to compare the computed exact to the approximate EIM.

Under the VGLM/VGAM approach, parameters can be modelled in terms of covariates. Particularly, if the standard deviation of the white noise is modelled in this way, then type.EIM = "exact" may certainly lead to unstable results. The reason is that white noise is a stationary process, and consequently, its variance must remain as a constant. Consequently, the use of variates to model this parameter contradicts the assumption of stationary random components to compute the exact EIMs proposed by Porat and Friedlander (1987).

To prevent convergence issues in such cases, this family function internally verifies whether the variance of the white noise remains as a constant at each Fisher scoring iteration. If this assumption is violated and type.EIM = "exact" is set, then AR1 automatically shifts to type.EIM = "approximate". Also, a warning is accordingly displayed.

Note

Multiple responses are handled. The mean is returned as the fitted values.

Author(s)

Victor Miranda (exact method) and Thomas W. Yee (approximate method).

References

Porat, B. and Friedlander, B. (1987). The Exact Cramer-Rao Bond for Gaussian Autoregressive Processes. IEEE Transactions on Aerospace and Electronic Systems, AES-23(4), 537–542.

See Also

AR1EIM, vglm.control, dAR1, arima.sim.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
### Example 1: using  arima.sim() to generate a 0-mean stationary time series.
nn <- 500
tsdata <- data.frame(x2 =  runif(nn))
ar.coef.1 <- rhobitlink(-1.55, inverse = TRUE)  # Approx -0.65
ar.coef.2 <- rhobitlink( 1.0, inverse = TRUE)   # Approx  0.50
set.seed(1)
tsdata  <- transform(tsdata,
              index = 1:nn,
              TS1 = arima.sim(nn, model = list(ar = ar.coef.1),
                              sd = exp(1.5)),
              TS2 = arima.sim(nn, model = list(ar = ar.coef.2),
                              sd = exp(1.0 + 1.5 * x2)))

### An autoregressive intercept--only model.   ###
### Using the exact EIM, and "nodrift = TRUE"  ###
fit1a <- vglm(TS1 ~ 1, data = tsdata, trace = TRUE,
              AR1(var.arg = FALSE, nodrift = TRUE,
                  type.EIM = "exact",
                  print.EIM = FALSE),
              crit = "coefficients")
Coef(fit1a)
summary(fit1a)

## Not run: 
### Two responses. Here, the white noise standard deviation of TS2   ###
### is modelled in terms of 'x2'. Also, 'type.EIM = exact'.  ###
fit1b <- vglm(cbind(TS1, TS2) ~ x2,
              AR1(zero = NULL, nodrift = TRUE,
                  var.arg = FALSE,
                  type.EIM = "exact"),
              constraints = list("(Intercept)" = diag(4),
                                 "x2" = rbind(0, 0, 1, 0)),
              data = tsdata, trace = TRUE, crit = "coefficients")
coef(fit1b, matrix = TRUE)
summary(fit1b)

### Example 2: another stationary time series
nn     <- 500
my.rho <- rhobitlink(1.0, inverse = TRUE)
my.mu  <- 1.0
my.sd  <- exp(1)
tsdata  <- data.frame(index = 1:nn, TS3 = runif(nn))

set.seed(2)
for (ii in 2:nn)
  tsdata$TS3[ii] <- my.mu/(1 - my.rho) +
                    my.rho * tsdata$TS3[ii-1] + rnorm(1, sd = my.sd)
tsdata <- tsdata[-(1:ceiling(nn/5)), ]  # Remove the burn-in data:

### Fitting an AR(1). The exact EIMs are used.
fit2a <- vglm(TS3 ~ 1, AR1(type.likelihood = "exact",  # "conditional",
                                type.EIM = "exact"),
              data = tsdata, trace = TRUE, crit = "coefficients")

Coef(fit2a)
summary(fit2a)      # SEs are useful to know

Coef(fit2a)["rho"]    # Estimate of rho, for intercept-only models
my.rho                # The 'truth' (rho)
Coef(fit2a)["drift"]  # Estimate of drift, for intercept-only models
my.mu /(1 - my.rho)   # The 'truth' (drift)

## End(Not run)

Example output

Loading required package: stats4
Loading required package: splines
VGLM    linear loop  1 :  coefficients =  1.5159735, -1.5484112
VGLM    linear loop  2 :  coefficients =  1.5158696, -1.5484224
VGLM    linear loop  3 :  coefficients =  1.5158577, -1.5484279
VGLM    linear loop  4 :  coefficients =  1.5158564, -1.5484303
VGLM    linear loop  5 :  coefficients =  1.5158562, -1.5484313
VGLM    linear loop  6 :  coefficients =  1.5158562, -1.5484318
VGLM    linear loop  7 :  coefficients =  1.5158562, -1.5484320
VGLM    linear loop  8 :  coefficients =  1.5158562, -1.5484321
        sd        rho 
 4.5533180 -0.6493743 

Call:
vglm(formula = TS1 ~ 1, family = AR1(var.arg = FALSE, nodrift = TRUE, 
    type.EIM = "exact", print.EIM = FALSE), data = tsdata, trace = TRUE, 
    crit = "coefficients")

Pearson residuals:
                    Min      1Q     Median     3Q   Max
loglink(sd)     -0.7083 -0.5925 -0.3956152 0.1958 8.663
rhobitlink(rho) -3.0308 -0.2819 -0.0007714 0.2480 4.668

Coefficients: 
              Estimate Std. Error z value Pr(>|z|)    
(Intercept):1  1.51586    0.02977   50.92   <2e-16 ***
(Intercept):2 -1.54843    0.08854  -17.49   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Names of linear predictors: loglink(sd), rhobitlink(rho)

Log-likelihood: -1467.671 on 998 degrees of freedom

Number of Fisher scoring iterations: 8 

No Hauck-Donner effect found in any of the estimates

VGLM    linear loop  1 :  coefficients = 
 1.51597353, -1.54841116,  1.42303496,  0.80284289,  1.11057123
VGLM    linear loop  2 :  coefficients = 
 1.51586961, -1.54842244,  1.03469736,  0.85622982,  1.56392433
VGLM    linear loop  3 :  coefficients = 
 1.51585774, -1.54842786,  1.01018433,  0.92417337,  1.58713679
VGLM    linear loop  4 :  coefficients = 
 1.51585638, -1.54843027,  1.05148929,  0.86576583,  1.53638731
VGLM    linear loop  5 :  coefficients = 
 1.51585623, -1.54843132,  0.98964717,  0.91190737,  1.60743688
VGLM    linear loop  6 :  coefficients = 
 1.51585621, -1.54843177,  1.09990833,  0.87505352,  1.47678562
VGLM    linear loop  7 :  coefficients = 
 1.51585621, -1.54843197,  0.93417133,  0.89635125,  1.67223450
VGLM    linear loop  8 :  coefficients = 
 1.51585621, -1.54843206,  1.25427681,  0.89614757,  1.29388328
VGLM    linear loop  9 :  coefficients = 
 1.51585621, -1.54843209,  0.89887599,  0.87743477,  1.71507935
VGLM    linear loop  10 :  coefficients = 
 1.51585621, -1.54843211,  1.38069688,  0.93126751,  1.14510620
VGLM    linear loop  11 :  coefficients = 
 1.51585621, -1.54843212,  0.97994001,  0.88500011,  1.62025697
VGLM    linear loop  12 :  coefficients = 
 1.51585621, -1.54843212,  1.12366705,  0.90461278,  1.44883994
VGLM    linear loop  13 :  coefficients = 
 1.51585621, -1.54843212,  0.91371349,  0.87793551,  1.69722926
VGLM    linear loop  14 :  coefficients = 
 1.51585621, -1.54843212,  1.32551893,  0.92692224,  1.20993065
VGLM    linear loop  15 :  coefficients = 
 1.51585621, -1.54843212,  0.93583119,  0.88063785,  1.67218248
VGLM    linear loop  16 :  coefficients = 
 1.51585621, -1.54843212,  1.24773817,  0.91799245,  1.30221503
VGLM    linear loop  17 :  coefficients = 
 1.51585621, -1.54843212,  0.89672080,  0.87636193,  1.71832383
VGLM    linear loop  18 :  coefficients = 
 1.51585621, -1.54843212,  1.38769457,  0.93358805,  1.13720096
VGLM    linear loop  19 :  coefficients = 
 1.51585621, -1.54843212,  0.98653503,  0.88623311,  1.61270117
VGLM    linear loop  20 :  coefficients = 
 1.51585621, -1.54843212,  1.10755400,  0.90224523,  1.46803040
VGLM    linear loop  21 :  coefficients = 
 1.51585621, -1.54843212,  0.92538662,  0.87947660,  1.68335036
VGLM    linear loop  22 :  coefficients = 
 1.5158562, -1.5484321,  1.2845460,  0.9219137,  1.2581966
VGLM    linear loop  23 :  coefficients = 
 1.51585621, -1.54843212,  0.91132284,  0.87787054,  1.70093732
VGLM    linear loop  24 :  coefficients = 
 1.51585621, -1.54843212,  1.33225459,  0.92753117,  1.20244815
VGLM    linear loop  25 :  coefficients = 
 1.51585621, -1.54843212,  0.94075766,  0.88121998,  1.66665726
VGLM    linear loop  26 :  coefficients = 
 1.51585621, -1.54843212,  1.23166121,  0.91618747,  1.32133961
VGLM    linear loop  27 :  coefficients = 
 1.51585621, -1.54843212,  0.89273415,  0.87594293,  1.72300659
VGLM    linear loop  28 :  coefficients = 
 1.5158562, -1.5484321,  1.4037534,  0.9353389,  1.1182747
VGLM    linear loop  29 :  coefficients = 
 1.51585621, -1.54843212,  1.00174129,  0.88788787,  1.59472269
VGLM    linear loop  30 :  coefficients = 
 1.5158562, -1.5484321,  1.0738955,  0.8984360,  1.5078536
Warning message:
In vglm.fitter(x = x, y = y, w = w, offset = offset, Xm2 = Xm2,  :
  convergence not obtained in 30 IRLS iterations
            loglink(sd1) rhobitlink(rho1) loglink(sd2) rhobitlink(rho2)
(Intercept)     1.515856        -1.548432     1.073895         0.898436
x2              0.000000         0.000000     1.507854         0.000000

Call:
vglm(formula = cbind(TS1, TS2) ~ x2, family = AR1(zero = NULL, 
    nodrift = TRUE, var.arg = FALSE, type.EIM = "exact"), data = tsdata, 
    constraints = list(`(Intercept)` = diag(4), x2 = rbind(0, 
        0, 1, 0)), trace = TRUE, crit = "coefficients")

Pearson residuals:
                     Min      1Q     Median     3Q   Max
loglink(sd1)     -0.7083 -0.5925 -0.3956152 0.1958 8.663
rhobitlink(rho1) -3.0308 -0.2819 -0.0007714 0.2480 4.668
loglink(sd2)     -1.9416 -0.7473 -0.3885596 0.1677 7.243
rhobitlink(rho2) -7.8475 -0.5002 -0.0074566 0.3791 9.923

Coefficients: 
              Estimate Std. Error z value Pr(>|z|)    
(Intercept):1  1.51586    0.02977  50.916   <2e-16 ***
(Intercept):2 -1.54843    0.08854 -17.489   <2e-16 ***
(Intercept):3  1.07390    0.10312  10.414   <2e-16 ***
(Intercept):4  0.89844    0.10242   8.772   <2e-16 ***
x2             1.50785    0.13574  11.109   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Names of linear predictors: loglink(sd1), rhobitlink(rho1), loglink(sd2), 
rhobitlink(rho2)

Log-likelihood: -3077.033 on 1995 degrees of freedom

Number of Fisher scoring iterations: 30 

No Hauck-Donner effect found in any of the estimates

VGLM    linear loop  1 :  coefficients = 
1.92023357, 0.99696614, 1.11291563
VGLM    linear loop  2 :  coefficients = 
1.92025423, 0.99686282, 1.11291324
VGLM    linear loop  3 :  coefficients = 
1.92025708, 0.99685211, 1.11291335
VGLM    linear loop  4 :  coefficients = 
1.92025680, 0.99685101, 1.11291335
VGLM    linear loop  5 :  coefficients = 
1.9202568, 0.9968509, 1.1129134
VGLM    linear loop  6 :  coefficients = 
1.92025681, 0.99685089, 1.11291335
    drift        sd       rho 
1.9202568 2.7097351 0.5053437 

Call:
vglm(formula = TS3 ~ 1, family = AR1(type.likelihood = "exact", 
    type.EIM = "exact"), data = tsdata, trace = TRUE, crit = "coefficients")

Pearson residuals:
                    Min      1Q    Median     3Q   Max
drift           -2.8733 -0.6789 -0.014500 0.7283 2.903
loglink(sd)     -0.6841 -0.5971 -0.339752 0.2148 4.777
rhobitlink(rho) -4.7050 -0.3711  0.005558 0.3789 3.104

Coefficients: 
              Estimate Std. Error z value Pr(>|z|)    
(Intercept):1  1.92026    0.13807  13.908   <2e-16 ***
(Intercept):2  0.99685    0.03352  29.736   <2e-16 ***
(Intercept):3  1.11291    0.11791   9.439   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Names of linear predictors: drift, loglink(sd), rhobitlink(rho)

Log-likelihood: -964.1564 on 1197 degrees of freedom

Number of Fisher scoring iterations: 6 

No Hauck-Donner effect found in any of the estimates

      rho 
0.5053437 
[1] 0.4621172
   drift 
1.920257 
[1] 1.859141

VGAM documentation built on Jan. 16, 2021, 5:21 p.m.