MARSS One-Step-Ahead Residuals

Calculates the standardized (or auxiliary) one-step-ahead residuals, aka the innovations residuals and their variance. Not exported. Access this function with MARSSresiduals(object, type="tt1"). To get the residuals as a data frame in long-form, use residuals(object, type="tt1").

Usage

MARSSresiduals.tt1(object, method = c("SS"), normalize = FALSE, 
    silent = FALSE, fun.kf = c("MARSSkfas", "MARSSkfss"))

Arguments

object: An object of class marssMLE.
method: Algorithm to use. Currently only "SS".
normalize: TRUE/FALSE See details.
silent: If TRUE, don't print inversion warnings.
fun.kf: Can be ignored. This will change the Kalman filter/smoother function from the value in object$fun.kf if desired.

Value

A list with the following components

model.residuals: The the observed one-step-ahead model residuals: data minus the model predictions conditioned on the data $t=1$ to $t-1$. These are termed innovations. A n x T matrix. NAs will appear where the data are missing.
state.residuals: The one-step-ahead state residuals $ \mathbf{x}_{t+1}^{t+1} - \mathbf{B}\mathbf{x}_{t}^t - \mathbf{u} $ . Note, state residual at time $t$ is the transition from time $t=t$ to $t+1$.
residuals: The residuals conditioned on the observed data up to time $t-1$. Returned as a (n+m) x T matrix with model.residuals in rows 1 to n and state.residuals in rows n+1 to n+m. NAs will appear in rows 1 to n in the places where data are missing.
var.residuals: The joint variance of the one-step-ahead residuals. Returned as a n+m x n+m x T matrix.
std.residuals: The Cholesky standardized residuals as a n+m x T matrix. This is residuals multiplied by the inverse of the lower triangle of the Cholesky decomposition of var.residuals. The model standardized residuals associated with the missing data are replaced with NA.
mar.residuals: The marginal standardized residuals as a n+m x T matrix. This is residuals multiplied by the inverse of the diagonal matrix formed by the square-root of the diagonal of var.residuals. The model marginal residuals associated with the missing data are replaced with NA.
bchol.residuals: The Block Cholesky standardized residuals as a (n+m) x T matrix. This is model.residuals multiplied by the inverse of the lower triangle of the Cholesky decomposition of var.residuals[1:n,1:n,] and state.residuals multiplied by the inverse of the lower triangle of the Cholesky decomposition of var.residuals[(n+1):(n+m),(n+1):(n+m),].
E.obs.residuals: The expected value of the model residuals conditioned on the observed data $t=1$ to $t-1$. Returned as a n x T matrix. Because all the data at time $t$ are unobserved for the purpose of estimation (since conditioning is from $t=1$ to $t-1$), this will be all 0s (unlike the case where we condition on the data from $t=1$ to $T$ or to $t$). This and var.obs.residuals are included for completeness since they are returned for MARSSresiduals.tT(), but they are not relevant for one-step-ahead residuals. See the discussion there.
var.obs.residuals: For one-step-ahead residuals, this will be the same as the 1:n, 1:n upper diagonal block in var.residuals since none of the $t$ data affect the residuals at time $t$ (the model residuals are conditioned only on the data up to $t-1$). This is different for smoothation residuals which are conditioned on the data from $t=1$ to $T$. This and E.obs.residuals are included for completeness since they are returned for MARSSresiduals.tT(), but they are not relevant for one-step-ahead residuals. See the discussion there. Note, also included as a code check. They are computed differently, but var.obs.residuals and var.residuals should always be the same.
msg: Any warning messages. This will be printed unless object$control$trace = -1 (suppress all error messages).

Details

This function returns the conditional expected value (mean) and variance of the one-step-ahead residuals. 'conditional' means in this context, conditioned on the observed data up to time $t-1$ and a set of parameters.

Model residuals

$\mathbf{v}_t$ is the difference between the data and the predicted data at time $t$ given $\mathbf{x}_t$: $$ \mathbf{v}_t = \mathbf{y}_t - \mathbf{Z} \mathbf{x}_t - \mathbf{a} - \mathbf{D}\mathbf{d}_t$$ The observed model residuals $\hat{\mathbf{v}}_t$ are the difference between the observed data and the predicted data at time $t$ using the fitted model. MARSSresiduals.tt1 fits the model using the data up to time $t-1$. So $$ \hat{\mathbf{v}}_t = \mathbf{y}_t - \mathbf{Z}\mathbf{x}_t^{t-1} - \mathbf{a} - \mathbf{D}\mathbf{d}_t$$ where $\mathbf{x}_t^{t-1}$ is the expected value of $\mathbf{X}_t$ conditioned on the data from $t=1$ to $t-1$ from the Kalman filter. $\mathbf{y}_t$ are your data and missing values will appear as NA.

State residuals

$\mathbf{w}_{t+1}$ are the difference between the state at time $t+1$ and the expected value of the state at time $t+1$ given the state at time $t$: $$ \mathbf{w}_{t+1} = \mathbf{x}_{t+1} - \mathbf{B} \mathbf{x}_{t} - \mathbf{u} - \mathbf{C}\mathbf{c}_{t+1}$$ The estimated state residuals $\hat{\mathbf{w}}_{t+1}$ are the difference between estimate of $\mathbf{x}_{t+1}$ minus the estimate using $\mathbf{x}_{t}$. $$ \hat{\mathbf{w}}_{t+1} = \mathbf{x}_{t+1}^{t+1} - \mathbf{B}\mathbf{x}_{t}^t - \mathbf{u} - \mathbf{C}\mathbf{c}_{t+1}$$ where $\mathbf{x}_{t+1}^{t+1}$ is the Kalman filter estimate of the states at time $t+1$ conditioned on the data up to time $t+1$ and $\mathbf{x}_{t}^t$ is the Kalman filter estimate of the states at time $t$ conditioned on the data up to time $t$. The estimated state residuals $\mathbf{w}_{t+1}$ are returned in state.residuals and rows $n+1$ to $n+m$ of residuals. state.residuals[,t] is $\mathbf{w}_{t+1}$ (notice time subscript difference). There are no NAs in the estimated state residuals (except for the last time step) as an estimate of the state exists whether or not there are associated data.

res1 and res2 in the code below will be the same.

dat <- t(harborSeal)[2:3,]
TT <- ncol(dat)
fit <- MARSS(dat)
B <- coef(fit, type="matrix")$B
U <- coef(fit, type="matrix")$U
xt <- MARSSkfss(fit)$xtt[,1:(TT-1)] # t 1 to TT-1
xtp1 <- MARSSkfss(fit)$xtt[,2:TT] # t 2 to TT
res1 <- xtp1 - B %*% xt - U %*% matrix(1,1,TT-1)
res2 <- MARSSresiduals(fit, type="tt1")$state.residuals

Joint residual variance

In a state-space model, $\mathbf{X}$ and $\mathbf{Y}$ are stochastic, and the model and state residuals are random variables $\hat{\mathbf{V}}_t$ and $\hat{\mathbf{W}}_{t+1}$. The joint distribution of $\hat{\mathbf{V}}_{t}, \hat{\mathbf{W}}_{t+1}$ is the distribution across all the different possible data sets that our MARSS equations with parameters $\Theta$ might generate. Denote the matrix of $\hat{\mathbf{V}}_{t}, \hat{\mathbf{W}}_{t+1}$, as $\widehat{\mathcal{E}}_{t}$. That distribution has an expected value (mean) and variance: $$ \textrm{E}[\widehat{\mathcal{E}}_t] = 0; \textrm{var}[\widehat{\mathcal{E}}_t] = \hat{\Sigma}_t $$ Our observed residuals residuals are one sample from this distribution. To standardize the observed residuals, we will use $ \hat{\Sigma}_t $. $ \hat{\Sigma}_t $ is returned in var.residuals. Rows/columns 1 to $n$ are the conditional variances of the model residuals and rows/columns $n+1$ to $n+m$ are the conditional variances of the state residuals. The off-diagonal blocks are the covariances between the two types of residuals. For one-step-ahead residuals (unlike smoothation residuals MARSSresiduals.tT), the covariance is zero.

var.residuals returned by this function is the conditional variance of the residuals conditioned on the data up to $t-1$ and the parameter set $\Theta$. The conditional variance for the model residuals is $$ \hat{\Sigma}_t = \mathbf{R}+\mathbf{Z}_t \mathbf{V}_t^{t-1} \mathbf{Z}_t^\top $$ where $\mathbf{V}_t^{t-1}$ is the variance of $\mathbf{X}_t$ conditioned on the data up to time $t-1$. This is returned by MARSSkf in Vtt1. The innovations variance is also returned in Sigma from MARSSkf and are used in the innovations form of the likelihood calculation.

Standardized residuals

std.residuals are Cholesky standardized residuals. These are the residuals multiplied by the inverse of the lower triangle of the Cholesky decomposition of the variance matrix of the residuals: $$ \hat{\Sigma}_t^{-1/2} \hat{\mathbf{v}}_t$$ These residuals are uncorrelated unlike marginal residuals.

The interpretation of the Cholesky standardized residuals is not straight-forward when the $\mathbf{Q}$ and $\mathbf{R}$ variance-covariance matrices are non-diagonal. The residuals which were generated by a non-diagonal variance-covariance matrices are transformed into orthogonal residuals in $\textrm{MVN}(0,\mathbf{I})$ space. For example, if v is 2x2 correlated errors with variance-covariance matrix R. The transformed residuals (from this function) for the i-th row of v is a combination of the row 1 effect and the row 1 effect plus the row 2 effect. So in this case, row 2 of the transformed residuals would not be regarded as solely the row 2 residual but rather how different row 2 is from row 1, relative to expected. If the errors are highly correlated, then the Cholesky standardized residuals can look rather non-intuitive.

mar.residuals are the marginal standardized residuals. These are the residuals multiplied by the inverse of the diagonal matrix formed from the square-root of the diagonal of the variance matrix of the residuals: $$ \textrm{dg}(\hat{\Sigma}_t)^{-1/2} \hat{\mathbf{v}}_t$$, where 'dg(A)' is the square matrix formed from the diagonal of A, aka diag(diag(A)). These residuals will be correlated if the variance matrix is non-diagonal.

The Block Cholesky standardized residuals are like the Cholesky standardized residuals except that the full variance-covariance matrix is not used, only the variance-covariance matrix for the model or state residuals (respectively) is used for standardization. For the one-step-ahead case, the model and state residuals are independent (unlike in the smoothations case) thus the Cholesky and Block Cholesky standardized residuals will be identical (unlike in the smoothations case).

Normalized residuals

If normalize=FALSE, the unconditional variance of $\mathbf{V}_t$ and $\mathbf{W}_t$ are $\mathbf{R}$ and $\mathbf{Q}$ and the model is assumed to be written as $$\mathbf{y}_t = \mathbf{Z} \mathbf{x}_t + \mathbf{a} + \mathbf{v}_t$$ $$\mathbf{x}_t = \mathbf{B} \mathbf{x}_{t-1} + \mathbf{u} + \mathbf{w}_t$$ If normalize=TRUE, the model is assumed to be written $$\mathbf{y}_t = \mathbf{Z} \mathbf{x}_t + \mathbf{a} + \mathbf{H}\mathbf{v}_t$$ $$\mathbf{x}_t = \mathbf{B} \mathbf{x}_{t-1} + \mathbf{u} + \mathbf{G}\mathbf{w}_t$$ with the variance of $\mathbf{V}_t$ and $\mathbf{W}_t$ equal to $\mathbf{I}$ (identity).

MARSSresiduals returns the residuals defined as in the first equations. To get the residuals defined as Harvey et al. (1998) define them (second equations), then use normalize=TRUE. In that case the unconditional variance of residuals will be $\mathbf{I}$ instead of $\mathbf{Q}$ and $\mathbf{R}$. Note, that the normalized residuals are not the same as the standardized residuals. In former, the unconditional residuals have a variance of $\mathbf{I}$ while in the latter it is the conditional residuals that have a variance of $\mathbf{I}$.

Author

Eli Holmes, NOAA, Seattle, USA.

References

R. H. Shumway and D. S. Stoffer (2006). Section on the calculation of the likelihood of state-space models in Time series analysis and its applications. Springer-Verlag, New York.

Holmes, E. E. 2014. Computation of standardized residuals for (MARSS) models. Technical Report. arXiv:1411.0045.

Examples

  dat <- t(harborSeal)
  dat <- dat[c(2,11),]
  fit <- MARSS(dat)
#> Success! abstol and log-log tests passed at 26 iterations.
#> Alert: conv.test.slope.tol is 0.5.
#> Test with smaller values (<0.1) to ensure convergence.
#> 
#> MARSS fit is
#> Estimation method: kem 
#> Convergence test: conv.test.slope.tol = 0.5, abstol = 0.001
#> Estimation converged in 26 iterations. 
#> Log-likelihood: 11.74016 
#> AIC: -9.480311   AICc: -6.3692   
#>  
#>                                           Estimate
#> R.diag                                      0.0115
#> U.X.CoastalEstuaries                        0.0613
#> U.X.OR.NorthCoast                           0.0510
#> Q.(X.CoastalEstuaries,X.CoastalEstuaries)   0.0147
#> Q.(X.OR.NorthCoast,X.OR.NorthCoast)         0.0122
#> x0.X.CoastalEstuaries                       7.3823
#> x0.X.OR.NorthCoast                          6.2707
#> Initial states (x0) defined at t=0
#> 
#> Standard errors have not been calculated. 
#> Use MARSSparamCIs to compute CIs and bias estimates.
#> 
  
  MARSSresiduals(fit, type="tt1")$std.residuals
#>                           [,1]         [,2]         [,3]     [,4]      [,5]
#> CoastalEstuaries   -0.05418676 -0.206286276  0.566365934 1.005768        NA
#> OR.NorthCoast               NA           NA -0.001882139       NA        NA
#> X.CoastalEstuaries -0.20628628  0.566365934  1.005767860 0.000000 0.2204673
#> X.OR.NorthCoast     0.00000000 -0.001882139  0.000000000 0.000000 0.0000000
#>                         [,6]       [,7]       [,8]       [,9]      [,10]
#> CoastalEstuaries   0.2204673  2.0405379  1.2425972 -0.7982010 -0.8746660
#> OR.NorthCoast             NA         NA -0.1395806  1.1554718  0.1779561
#> X.CoastalEstuaries 2.0405379  1.2425972 -0.7982010 -0.8746660  1.3133883
#> X.OR.NorthCoast    0.0000000 -0.1395806  1.1554718  0.1779561  0.3438440
#>                         [,11]      [,12]      [,13]      [,14]       [,15]
#> CoastalEstuaries    1.3133883 -1.1010233  1.7849809 -0.3351368 -0.71142906
#> OR.NorthCoast       0.3438440 -0.8332411  1.6166181  0.6909545  1.09517141
#> X.CoastalEstuaries -1.1010233  1.7849809 -0.3351368 -0.7114291  0.00000000
#> X.OR.NorthCoast    -0.8332411  1.6166181  0.6909545  1.0951714  0.03597429
#>                          [,16]      [,17]       [,18]       [,19]      [,20]
#> CoastalEstuaries            NA  0.5726369 -0.76775282 -0.31659312 -2.3132481
#> OR.NorthCoast       0.03597429 -0.9304752  1.64162505 -0.08266883 -0.6865966
#> X.CoastalEstuaries  0.57263695 -0.7677528 -0.31659312 -2.31324808 -0.4943444
#> X.OR.NorthCoast    -0.93047525  1.6416251 -0.08266883 -0.68659660 -1.3199941
#>                          [,21]       [,22]      [,23]      [,24]      [,25]
#> CoastalEstuaries   -0.49434441  0.05447392 0.01643296         NA -0.8768965
#> OR.NorthCoast      -1.31999410 -0.07484734         NA  0.2166126 -2.1429904
#> X.CoastalEstuaries  0.05447392  0.01643296 0.00000000 -0.8768965  0.0000000
#> X.OR.NorthCoast    -0.07484734  0.00000000 0.21661261 -2.1429904 -1.5204468
#>                         [,26]      [,27]     [,28]     [,29] [,30]
#> CoastalEstuaries           NA         NA        NA        NA    NA
#> OR.NorthCoast      -1.5204468 -0.8263151 0.9151561 0.7476442    NA
#> X.CoastalEstuaries  0.0000000  0.0000000 0.0000000 0.0000000    NA
#> X.OR.NorthCoast    -0.8263151  0.9151561 0.7476442 0.0000000    NA
  residuals(fit, type="tt1")
#>    type        .rownames  name  t    value  .fitted       .resids    .sigma
#> 1  ytt1 CoastalEstuaries model  1 7.434848 7.443613 -0.0087655032 0.1617647
#> 2  ytt1 CoastalEstuaries model  2 7.462789 7.500042 -0.0372528292 0.1805880
#> 3  ytt1 CoastalEstuaries model  3 7.641084 7.537255  0.1038290197 0.1833250
#> 4  ytt1 CoastalEstuaries model  4 7.851661 7.666950  0.1847108729 0.1836516
#> 5  ytt1 CoastalEstuaries model  5       NA 7.850112            NA        NA
#> 6  ytt1 CoastalEstuaries model  6 7.959975 7.911459  0.0485152957 0.2200567
#> 7  ytt1 CoastalEstuaries model  7 8.391176 8.009815  0.3813609196 0.1868923
#> 8  ytt1 CoastalEstuaries model  8 8.555837 8.327130  0.2287064413 0.1840552
#> 9  ytt1 CoastalEstuaries model  9 8.392990 8.539648 -0.1466584134 0.1837362
#> 10 ytt1 CoastalEstuaries model 10 8.343554 8.504229 -0.1606755891 0.1836994
#> 11 ytt1 CoastalEstuaries model 11 8.700847 8.459584  0.2412630005 0.1836951
#> 12 ytt1 CoastalEstuaries model 12 8.477828 8.680081 -0.2022520356 0.1836946
#> 13 ytt1 CoastalEstuaries model 13 8.935904 8.608012  0.3278912699 0.1836945
#> 14 ytt1 CoastalEstuaries model 14 8.824089 8.885652 -0.0615627982 0.1836945
#> 15 ytt1 CoastalEstuaries model 15 8.775704 8.906390 -0.1306856341 0.1836945
#> 16 ytt1 CoastalEstuaries model 16       NA 8.881530            NA        NA
#> 17 ytt1 CoastalEstuaries model 17 9.068892 8.942877  0.1260149752 0.2200609
#> 18 ytt1 CoastalEstuaries model 18 8.956866 9.100354 -0.1434873391 0.1868926
#> 19 ytt1 CoastalEstuaries model 19 9.007122 9.065393 -0.0582706082 0.1840552
#> 20 ytt1 CoastalEstuaries model 20 8.663196 9.088224 -0.4250273931 0.1837362
#> 21 ytt1 CoastalEstuaries model 21 8.778326 8.869137 -0.0908107534 0.1836994
#> 22 ytt1 CoastalEstuaries model 22 8.880586 8.870579  0.0100065927 0.1836951
#> 23 ytt1 CoastalEstuaries model 23 8.941545 8.938527  0.0030186458 0.1836946
#> 24 ytt1 CoastalEstuaries model 24       NA 9.001865            NA        NA
#> 25 ytt1 CoastalEstuaries model 25 8.870242 9.063212 -0.1929705976 0.2200609
#> 26 ytt1 CoastalEstuaries model 26       NA 8.977353            NA        NA
#> 27 ytt1 CoastalEstuaries model 27       NA 9.038700            NA        NA
#> 28 ytt1 CoastalEstuaries model 28       NA 9.100047            NA        NA
#> 29 ytt1 CoastalEstuaries model 29       NA 9.161394            NA        NA
#> 30 ytt1 CoastalEstuaries model 30       NA 9.222741            NA        NA
#> 31 ytt1    OR.NorthCoast model  1       NA 6.321668            NA        NA
#> 32 ytt1    OR.NorthCoast model  2       NA 6.372664            NA        NA
#> 33 ytt1    OR.NorthCoast model  3 6.423247 6.423659 -0.0004124727 0.2191510
#> 34 ytt1    OR.NorthCoast model  4       NA 6.474341            NA        NA
#> 35 ytt1    OR.NorthCoast model  5       NA 6.525337            NA        NA
#> 36 ytt1    OR.NorthCoast model  6       NA 6.576333            NA        NA
#> 37 ytt1    OR.NorthCoast model  7       NA 6.627329            NA        NA
#> 38 ytt1    OR.NorthCoast model  8 6.638568 6.678324 -0.0397565172 0.2848283
#> 39 ytt1    OR.NorthCoast model  9 6.906755 6.695192  0.2115631430 0.1830968
#> 40 ytt1    OR.NorthCoast model 10 6.916715 6.885274  0.0314412894 0.1766801
#> 41 ytt1    OR.NorthCoast model 11 7.016610 6.956143  0.0604665822 0.1758547
#> 42 ytt1    OR.NorthCoast model 12 6.898715 7.045150 -0.1464351316 0.1757416
#> 43 ytt1    OR.NorthCoast model 13 7.288244 7.004163  0.2840818171 0.1757260
#> 44 ytt1    OR.NorthCoast model 14 7.355002 7.233585  0.1214171679 0.1757238
#> 45 ytt1    OR.NorthCoast model 15 7.553287 7.360839  0.1924473895 0.1757235
#> 46 ytt1    OR.NorthCoast model 16 7.539027 7.532706  0.0063215283 0.1757235
#> 47 ytt1    OR.NorthCoast model 17 7.424165 7.587672 -0.1635063539 0.1757235
#> 48 ytt1    OR.NorthCoast model 18 7.824446 7.535974  0.2884720744 0.1757235
#> 49 ytt1    OR.NorthCoast model 19 7.753624 7.768150 -0.0145268555 0.1757235
#> 50 ytt1    OR.NorthCoast model 20 7.689371 7.810022 -0.1206511466 0.1757235
#> 51 ytt1    OR.NorthCoast model 21 7.553287 7.785241 -0.2319539627 0.1757235
#> 52 ytt1    OR.NorthCoast model 22 7.677400 7.690553 -0.0131524352 0.1757235
#> 53 ytt1    OR.NorthCoast model 23       NA 7.733288            NA        NA
#> 54 ytt1    OR.NorthCoast model 24 7.829233 7.784284  0.0449488452 0.2075080
#> 55 ytt1    OR.NorthCoast model 25 7.484369 7.868240 -0.3838710165 0.1791287
#> 56 ytt1    OR.NorthCoast model 26 7.404888 7.672761 -0.2678732130 0.1761806
#> 57 ytt1    OR.NorthCoast model 27 7.409742 7.554997 -0.1452549981 0.1757865
#> 58 ytt1    OR.NorthCoast model 28 7.675546 7.514724  0.1608223910 0.1757322
#> 59 ytt1    OR.NorthCoast model 29 7.798113 7.666733  0.1313795472 0.1757247
#> 60 ytt1    OR.NorthCoast model 30       NA 7.800245            NA        NA
#>     .std.resids
#> 1  -0.054186758
#> 2  -0.206286276
#> 3   0.566365934
#> 4   1.005767860
#> 5            NA
#> 6   0.220467252
#> 7   2.040537901
#> 8   1.242597244
#> 9  -0.798201023
#> 10 -0.874666007
#> 11  1.313388330
#> 12 -1.101023270
#> 13  1.784980949
#> 14 -0.335136785
#> 15 -0.711429056
#> 16           NA
#> 17  0.572636946
#> 18 -0.767752824
#> 19 -0.316593118
#> 20 -2.313248083
#> 21 -0.494344408
#> 22  0.054473923
#> 23  0.016432959
#> 24           NA
#> 25 -0.876896499
#> 26           NA
#> 27           NA
#> 28           NA
#> 29           NA
#> 30           NA
#> 31           NA
#> 32           NA
#> 33 -0.001882139
#> 34           NA
#> 35           NA
#> 36           NA
#> 37           NA
#> 38 -0.139580626
#> 39  1.155471842
#> 40  0.177956078
#> 41  0.343844019
#> 42 -0.833241083
#> 43  1.616618091
#> 44  0.690954477
#> 45  1.095171415
#> 46  0.035974293
#> 47 -0.930475250
#> 48  1.641625053
#> 49 -0.082668833
#> 50 -0.686596598
#> 51 -1.319994101
#> 52 -0.074847339
#> 53           NA
#> 54  0.216612615
#> 55 -2.142990433
#> 56 -1.520446816
#> 57 -0.826315058
#> 58  0.915156126
#> 59  0.747644217
#> 60           NA