Residuals of a statespace system

This internal helper function computes the residuals (and the directional derivatives of the residuals) for a statespace system of the form $$a_{t+1} = A a_t + B u_t, \; y_t = C a_t + D u_t$$ The system must be square and non-empty, i.e. \(m=n>0\).

Usage

residuals_STSP_cpp(A, B, C, D, y, a, u, dPI, dU)

Arguments

A

\((s,s)\) matrix.

B

\((s,m)\) matrix.

C

\((m,s)\) matrix.

D

\((m,m)\) matrix, must be regular.

y

\((m,N)\) matrix with the outputs: \((y_1,y_2,\ldots,y_N)\).

a

\((s,N+1)\) matrix. This matrix is overwritten with the (computed) states: \((a_1,a_2,\ldots,a_N,a_{N+1})\). On input a[,1] must hold the initial state \(a_1\).

u

\((m,N)\) matrix. This matrix is overwritten with (computed) residuals: \((u_1,u_2,\ldots,u_N)\).

dPI

\(((m+s)^2,K)\) matrix.

dU

\((mN,K)\) matrix or \((0,0)\) matrix. This matrix is overwritten with the directional derivatives of the residuals. However, if the matrix is empty then no derivatives are computed.

Value

This RcppArmadillo implementation returns NULL but overwrites the input arguments a, u and dU!

Note

Use this procedure with care!

• The procedure does not check the input arguments.

• The procedure overwrites the input arguments a, u and dU.

• The data matrices are organized columnwise (to avoid memory shuffling)!

See also

outputs_ARMA_cpp, residuals_ARMA_cpp, cll_theta_ARMA_cpp, outputs_STSP_cpp, residuals_STSP_cpp, cll_theta_STSP_cpp and solve_de, solve_inverse_de and ll.

Examples

# generate a random statespace model (3 outputs, 3 inputs and 4 states)
m = 2
s = 3
model = test_stspmod(dim = c(m, m), s = s, digits = 2)

# generate random noise sequence (sample size N = 10)
n.obs = 10
u = matrix(rnorm(n.obs*m), nrow = m, ncol = n.obs)

# generate matrix for the state sequence
a = matrix(0, nrow = s, ncol = n.obs+1)
a[,1] = rnorm(s) # random initial state a[1]

# generate matrix for the outputs
y = matrix(0, nrow = m, ncol = n.obs)

# compute outputs and states
outputs_STSP_cpp(model$sys$A, model$sys$B, model$sys$C, model$sys$D, u, a, y)

# recompute the states and disturbances/residuals from the given outputs:
uu = u + 0 # "deep copy" of the disturbances
aa = a + 0 # and the states
aa[, 2:(n.obs+1)] = 0 # clear all states a[t], t > 1
residuals_STSP_cpp(model$sys$A, model$sys$B, model$sys$C, model$sys$D,
                   y, aa, uu, diag(0), diag(0))
all.equal(u, uu) # check
#> [1] TRUE
all.equal(a, aa) # check
#> [1] TRUE

# compute directional derivatives of residuals
dPI = diag((m+s)^2)
dU = matrix(0, nrow = n.obs*m, ncol = ncol(dPI))
residuals_STSP_cpp(model$sys$A, model$sys$B, model$sys$C, model$sys$D,
                   y, aa, uu, dPI, dU)

# check the directional derivatives
eps = 1e-8
dU_num = matrix(0, nrow = m*n.obs, ncol = (m+s)^2)
dPI = matrix(0, nrow = (m+s), ncol = (m+s))
for (k in (1:((m+s)^2))) {
  dPI[] = 0
  dPI[k] = eps
  uu = u + 0 # "deep copy" of the disturbances
  aa = a + 0 # and the states
  residuals_STSP_cpp(model$sys$A + dPI[1:s,1:s],
                     model$sys$B + dPI[1:s,(s+1):(s+m)],
                     model$sys$C + dPI[(s+1):(s+m),1:s],
                     model$sys$D + dPI[(s+1):(s+m),(s+1):(s+m)],
                     y, aa, uu, diag(0), diag(0))
  dU_num[, k] = c(uu - u )/eps  # num. approx. of the derivative in direction "dPI"
}
# relative error of the numerical approximation
junk = (abs(dU)+abs(dU_num))
junk[junk == 0] = 1
2*abs(dU_num - dU)/junk
#>               [,1]         [,2]         [,3]         [,4]         [,5]
#>  [1,] 2.000000e+00 2.000000e+00 2.000000e+00 1.075932e-08 2.000000e+00
#>  [2,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 1.075932e-08
#>  [3,] 3.186212e-08 3.579097e-08 4.039238e-08 1.752661e-08 5.803432e-08
#>  [4,] 1.074962e-07 1.489792e-07 3.852173e-08 3.460251e-08 1.598594e-07
#>  [5,] 1.157524e-07 5.095788e-08 2.537305e-07 5.217776e-09 3.326926e-08
#>  [6,] 1.315437e-07 3.374915e-07 5.377632e-08 6.339870e-07 4.053750e-08
#>  [7,] 6.685966e-08 6.484088e-07 6.607462e-08 6.397383e-08 5.599988e-08
#>  [8,] 1.664661e-07 2.324767e-06 8.802958e-09 1.160632e-07 6.577718e-09
#>  [9,] 3.479932e-07 1.844783e-07 1.162988e-08 5.336204e-08 7.423279e-09
#> [10,] 5.741504e-07 1.071168e-07 9.153993e-09 2.286567e-08 2.682838e-08
#> [11,] 5.049264e-08 6.576649e-08 8.336142e-09 1.792653e-08 4.070700e-08
#> [12,] 5.960548e-08 1.806476e-07 2.004137e-08 6.831834e-08 2.439828e-08
#> [13,] 3.995683e-08 1.183683e-08 2.742813e-08 1.405923e-08 2.137221e-08
#> [14,] 1.323205e-07 2.033844e-06 8.872500e-09 1.151428e-06 2.016214e-08
#> [15,] 7.170350e-08 2.601382e-07 1.185412e-08 8.475312e-08 2.838888e-08
#> [16,] 1.619748e-07 1.594132e-07 2.916319e-08 3.897970e-08 3.663698e-08
#> [17,] 1.665679e-08 5.132484e-07 9.552205e-09 3.685681e-08 2.991595e-08
#> [18,] 5.070742e-08 1.041959e-07 3.124539e-08 5.545476e-08 3.925251e-08
#> [19,] 5.095013e-08 2.702853e-08 4.421457e-08 4.771770e-09 4.753418e-08
#> [20,] 4.661526e-08 1.929875e-06 2.001892e-08 9.108789e-07 3.605221e-08
#>               [,6]         [,7]         [,8]         [,9]        [,10]
#>  [1,] 2.000000e+00 2.000000e+00 2.000000e+00 1.131678e-06 2.000000e+00
#>  [2,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 1.113355e-07
#>  [3,] 8.170210e-08 1.034502e-06 1.413455e-07 5.321547e-09 4.466609e-07
#>  [4,] 2.704707e-06 3.855674e-06 1.182311e-06 1.140510e-06 3.185600e-08
#>  [5,] 4.623740e-08 4.183661e-08 2.338541e-08 7.093364e-08 5.508330e-09
#>  [6,] 1.356550e-07 1.097172e-07 3.029195e-08 5.278920e-08 3.648631e-08
#>  [7,] 1.664981e-07 2.267173e-08 3.670966e-08 4.388042e-06 3.731196e-08
#>  [8,] 5.022520e-08 6.556526e-08 5.763200e-09 1.959923e-08 1.138926e-08
#>  [9,] 3.141253e-06 1.063464e-07 5.437144e-09 1.888483e-08 2.545790e-08
#> [10,] 1.928351e-08 2.601412e-08 7.670408e-10 4.794793e-08 1.792162e-09
#> [11,] 8.347454e-09 9.221605e-08 3.226239e-09 8.947903e-08 2.491142e-11
#> [12,] 2.763057e-07 7.455582e-08 1.153922e-08 3.294806e-08 4.381306e-09
#> [13,] 8.993611e-08 2.257600e-08 2.171724e-08 9.296685e-09 2.771293e-09
#> [14,] 2.052089e-07 6.968366e-08 1.350649e-09 3.507194e-08 8.978265e-09
#> [15,] 4.348307e-07 7.688941e-08 1.056695e-08 2.651263e-08 1.867296e-08
#> [16,] 4.675003e-08 1.189991e-07 1.295404e-08 7.205318e-08 2.107867e-09
#> [17,] 2.834989e-08 8.952275e-08 7.284231e-09 5.667003e-08 1.110488e-09
#> [18,] 9.027183e-08 1.856026e-07 1.196378e-08 8.227280e-08 1.248393e-09
#> [19,] 3.722311e-08 4.330361e-07 1.380998e-08 1.732444e-07 6.872835e-09
#> [20,] 2.404268e-07 6.303661e-08 3.218090e-09 3.446218e-08 1.898721e-10
#>              [,11]        [,12]        [,13]        [,14]        [,15]
#>  [1,] 2.000000e+00 2.000000e+00 2.000000e+00 1.086048e-09 2.000000e+00
#>  [2,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 2.513925e-08
#>  [3,] 3.315354e-10 1.889513e-08 9.570703e-09 1.454498e-07 2.977446e-08
#>  [4,] 6.182345e-08 1.455513e-07 6.177393e-10 2.954215e-08 2.392999e-08
#>  [5,] 7.145524e-08 5.885461e-08 1.330320e-07 1.205470e-07 7.714745e-08
#>  [6,] 7.668011e-08 1.559413e-07 2.364268e-08 3.020189e-10 3.290929e-08
#>  [7,] 3.067631e-08 7.092089e-08 1.313786e-07 6.034286e-09 1.380139e-07
#>  [8,] 4.410584e-08 1.006175e-06 2.056192e-08 2.225353e-07 1.448901e-08
#>  [9,] 6.070601e-08 1.288510e-06 9.924175e-10 1.855999e-06 1.081766e-08
#> [10,] 4.754040e-08 2.572919e-06 4.787631e-08 2.837652e-08 2.905774e-08
#> [11,] 8.540794e-09 9.529981e-08 2.680301e-08 1.631147e-08 7.087777e-08
#> [12,] 5.377599e-07 1.833774e-07 2.689266e-08 5.390497e-09 1.093626e-08
#> [13,] 2.898540e-09 9.058536e-08 2.656695e-08 1.374390e-08 2.866201e-08
#> [14,] 6.094147e-08 3.283779e-07 1.098062e-08 1.397911e-07 5.900570e-09
#> [15,] 4.505867e-08 2.314631e-07 1.900022e-08 6.981666e-08 1.466670e-08
#> [16,] 1.805968e-07 6.958001e-07 8.494064e-09 1.734512e-07 2.478770e-08
#> [17,] 7.105628e-08 9.023277e-08 1.412866e-08 2.368654e-08 3.578909e-08
#> [18,] 1.910255e-07 2.309797e-07 6.695830e-09 2.307009e-08 5.475946e-09
#> [19,] 1.579230e-08 4.566930e-07 4.311102e-09 2.320777e-09 4.454991e-08
#> [20,] 6.451751e-08 2.064840e-07 3.817754e-09 5.819332e-08 3.805869e-09
#>              [,16]        [,17]        [,18]        [,19]        [,20]
#>  [1,] 2.000000e+00 2.000000e+00 2.000000e+00 5.899136e-09 2.000000e+00
#>  [2,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 1.176494e-08
#>  [3,] 9.718396e-09 2.438214e-08 9.016326e-09 1.061842e-08 1.823861e-08
#>  [4,] 8.221098e-08 1.420344e-07 3.554955e-08 1.318116e-09 1.400125e-07
#>  [5,] 2.331672e-08 6.331857e-07 2.950538e-07 1.635135e-08 5.070854e-07
#>  [6,] 7.408486e-08 3.087400e-07 2.975985e-08 2.056495e-08 9.025849e-08
#>  [7,] 7.669562e-08 1.722099e-07 3.167598e-08 8.144718e-09 6.545394e-08
#>  [8,] 1.592914e-07 1.338047e-07 2.342117e-09 1.062973e-07 2.319386e-07
#>  [9,] 5.685784e-08 2.720696e-07 2.510647e-07 4.907111e-08 1.385100e-07
#> [10,] 2.151694e-08 6.982250e-07 4.593984e-07 1.218809e-07 6.646196e-07
#> [11,] 7.055892e-09 9.588819e-07 5.579498e-07 2.831668e-07 3.640883e-06
#> [12,] 3.107749e-08 4.201407e-08 4.885709e-09 4.784480e-08 2.418832e-07
#> [13,] 2.286579e-07 3.114189e-08 6.753980e-08 4.876711e-09 1.988561e-07
#> [14,] 1.519367e-07 6.213371e-07 1.277538e-07 2.430001e-07 2.876687e-07
#> [15,] 1.255598e-07 8.545111e-07 2.068658e-07 8.630924e-08 2.250087e-07
#> [16,] 2.440180e-07 1.959124e-06 2.965566e-08 1.082626e-06 7.790226e-07
#> [17,] 1.459087e-07 2.175077e-07 1.174725e-07 6.894526e-07 1.452562e-06
#> [18,] 3.304012e-07 2.127956e-07 8.240525e-08 1.927536e-07 4.441959e-07
#> [19,] 8.781336e-08 5.216887e-08 8.231572e-08 1.282847e-07 6.320619e-07
#> [20,] 2.029027e-07 1.551358e-06 1.607406e-07 7.666133e-07 3.202635e-07
#>              [,21]        [,22]        [,23]        [,24]        [,25]
#>  [1,] 2.000000e+00 2.000000e+00 2.000000e+00 6.516317e-09 2.000000e+00
#>  [2,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 2.071710e-08
#>  [3,] 6.245522e-08 1.689095e-08 7.901373e-09 1.809532e-08 3.038508e-08
#>  [4,] 1.371075e-07 2.689901e-07 5.543695e-09 5.011156e-09 2.006105e-08
#>  [5,] 9.563963e-08 1.682724e-06 6.888426e-07 5.426523e-09 1.016955e-06
#>  [6,] 1.142194e-07 4.584656e-07 6.200320e-08 1.060949e-07 2.168329e-08
#>  [7,] 6.464573e-08 6.537247e-07 3.637412e-08 2.683827e-07 6.871651e-08
#>  [8,] 4.547309e-08 1.192176e-06 7.989220e-08 4.254949e-07 5.568590e-08
#>  [9,] 2.031046e-07 2.550322e-07 2.272690e-07 3.731935e-08 1.075565e-07
#> [10,] 8.515991e-07 7.017187e-07 3.455242e-08 4.384980e-07 1.788479e-07
#> [11,] 1.651453e-07 7.061197e-06 1.519340e-07 4.827146e-07 5.140289e-07
#> [12,] 3.571436e-08 6.723812e-07 3.776869e-08 3.469672e-07 8.855761e-08
#> [13,] 1.522266e-07 3.470914e-07 6.861230e-08 2.444851e-07 1.210583e-07
#> [14,] 1.324700e-07 7.245459e-06 1.176127e-07 3.546585e-06 4.555955e-08
#> [15,] 1.035235e-07 3.363364e-04 2.900984e-07 3.542079e-07 9.460380e-08
#> [16,] 3.572550e-06 2.046287e-07 7.740472e-08 1.157556e-06 1.533776e-07
#> [17,] 5.842137e-07 3.004976e-06 1.099562e-07 2.796271e-06 3.178225e-07
#> [18,] 3.822025e-08 7.988323e-07 1.302892e-07 6.505714e-07 8.016194e-08
#> [19,] 4.667872e-07 2.231551e-07 7.938793e-08 5.427415e-07 1.673795e-07
#> [20,] 9.401545e-08 2.726302e-06 2.054776e-07 1.781893e-06 5.233476e-08