Ho-Kalman Realization Algorithm — pseries2stsp • rationalmatrices

This helper function implements the Ho-Kalman algorithm.

Usage

pseries2stsp(
  obj,
  method = c("balanced", "echelon"),
  Hsize = NULL,
  s = NULL,
  nu = NULL,
  tol = sqrt(.Machine$double.eps),
  Wrow = NULL,
  Wcol = NULL
)

Arguments

obj: pseries object or 3-D array with dimension \((m,n,l+1)\).
method: Character string, which determines the method and the "parametrization" type of the state space model. See below for more details.
Hsize: integer vector c(f,p), number of block rows and block columns of the Hankel matrix which is used to construct the statespace realization. If NULL a default choice is made.
s: desired state dimension. Only used for method = "balanced". Note however, if \(s\) is larger than the rank of the Hankel matrix, then the procedure will break down. If s is missing, then the state dimension is determined from the singular values of the Hankel matrix. To be precise the state dimension is chosen as the number of singular values which are greater than or equal to tol times the maximum singular value.
nu: Kronecker indices. Only used for method = "echelon". If missing, then nu is computed with a QR decomposition of the transpose of the Hankel matrix of the impulse response coefficients.
tol: tolerance parameter used for the QR decomposition or the SVD decomposition of the Hankel matrix \(H\) of the impulse response coefficients.
Wrow, Wcol: weighting matrices (default is no weighting, i.e. identity matrices). These weighting matrices are only used for method="balanced", where the SVD of the weighted Hankel matrix Wrow %*% H %*% t(Wcol) is computed.

Value

List with slots

Xs: stsp object, the rational matrix in statespace form
Hsv: Singular values of the Hankel matrix for method='balanced' and NULL else.
nu: Kronecker indices for method='echelon' and NULL else.

Details

The procedure(s) may be used for model reduction (with some care).

There are a number of restrictions on the number of lags \(l\) of the impulse response, the number of block rows (\(f\)), block columns (\(p\)) of the Hankel matrix and the Kronecker indices \(\nu_i\). We require that: \(p>0\), \(f>1\), \(l \geq f+p-1\) and \(\nu_i <f\). If these restrictions are not satisfied an error is thrown.

Examples

# generate random rational matrix X(z) in statespace form
# make sure that the A matrix is stable
m = 3
n = 2
s = 7
A = matrix(rnorm(s*s), nrow = s, ncol = s)
A = A / (1.1 * max(abs(eigen(A, only.values = TRUE)$values)))
Xs = stsp(A, B = matrix(rnorm(s*n), nrow = s, ncol = n),
          C = matrix(rnorm(s*m), nrow = m, ncol = s),
          D = diag(1, nrow = m, ncol = n))
Xi = pseries(Xs, lag.max = 20)

out = pseries2stsp(Xi, method = 'balanced')
print(out)
#> $Xs
#> statespace realization [3,2] with s = 7 states
#>               s[1]         s[2]         s[3]        s[4]         s[5]
#> s[1]  0.3313851522  0.770307985 -0.111460220 -0.09982590 -0.066785063
#> s[2] -0.7921919615  0.526516640 -0.034292547  0.07457127  0.102487584
#> s[3] -0.2078866452 -0.267664979 -0.389684672 -0.17581638 -0.078973283
#> s[4] -0.0805892165 -0.032096117  0.431976222 -0.33424645  0.183531989
#> s[5] -0.0154308669 -0.018511367  0.193517532  0.66435507 -0.329885130
#> s[6]  0.0001131816 -0.001228138  0.008490031  0.17651995  0.506428961
#> s[7]  0.0198225494 -0.001616940  0.044497523 -0.10905861 -0.008773808
#> x[1] -0.3151281853  0.529051907  0.347009515  0.39861346 -0.571473592
#> x[2] -2.0953837797 -0.684634712 -0.448603366 -0.42842455 -0.475306808
#> x[3]  0.5406479391 -0.258394568 -1.809499404  0.39147902  0.154531626
#>             s[6]          s[7]        u[1]        u[2]
#> s[1] -0.01877308  0.0001067318 -0.31496301 -2.23977296
#> s[2]  0.01172167  0.0095199948  0.81048633  0.56074243
#> s[3]  0.17315009  0.0569796953  0.30973252 -1.55337413
#> s[4] -0.01768040  0.1708131766 -0.42164705 -0.33498274
#> s[5] -0.14636226  0.1662406172  0.22956890 -0.21700855
#> s[6]  0.23419911 -0.2582851157  0.06007157 -0.04202095
#> s[7]  0.02504584 -0.2590689778  0.29512327 -0.01139882
#> x[1]  0.27998343 -0.0471096862  1.00000000  0.00000000
#> x[2] -0.24640435 -0.0703576446  0.00000000  1.00000000
#> x[3] -0.07471873  0.0248091034  0.00000000  0.00000000
#> 
#> $Hsv
#>  [1] 1.433126e+01 1.282354e+01 4.840859e+00 1.496822e+00 1.052230e+00
#>  [6] 3.378444e-01 1.283687e-01 1.449573e-15 1.147411e-15 8.610457e-16
#> [11] 6.310861e-16 4.615842e-16 4.306397e-16 3.479924e-16 3.183540e-16
#> [16] 2.397370e-16 2.024793e-16 1.417466e-16 1.266129e-16 1.167343e-16
#> 
#> $nu
#> NULL
#> 
# check impulse response
all.equal(pseries(out$Xs, lag.max = 20), Xi)
#> [1] TRUE

Xs1 = as.stsp(Xi)
all.equal(Xs1, out$Xs)
#> [1] TRUE

out = pseries2stsp(Xi, method = 'echelon')
print(out)
#> $Xs
#> statespace realization [3,2] with s = 7 states
#>            s[1]        s[2]       s[3]       s[4]       s[5]       s[6]
#> s[1]  0.0000000  0.00000000  0.0000000  1.0000000  0.0000000  0.0000000
#> s[2]  0.0000000  0.00000000  0.0000000  0.0000000  1.0000000  0.0000000
#> s[3]  0.0000000  0.00000000  0.0000000  0.0000000  0.0000000  1.0000000
#> s[4]  0.0000000  0.00000000  0.0000000  0.0000000  0.0000000  0.0000000
#> s[5]  0.2301243 -0.35973214  0.1484191 -1.0397365  1.0968456 -0.2250004
#> s[6] -1.0493434  0.06555365 -0.7206394 -1.0323758 -1.2648783 -1.3069378
#> s[7] -0.7433376 -0.20534077 -0.5408140 -0.6749443 -0.1149876 -0.6724519
#> x[1]  1.0000000  0.00000000  0.0000000  0.0000000  0.0000000  0.0000000
#> x[2]  0.0000000  1.00000000  0.0000000  0.0000000  0.0000000  0.0000000
#> x[3]  0.0000000  0.00000000  1.0000000  0.0000000  0.0000000  0.0000000
#>             s[7]         u[1]       u[2]
#> s[1]  0.00000000  0.339172322  0.4426997
#> s[2]  0.00000000  0.002096456  5.2639457
#> s[3]  0.00000000 -1.066927231  1.2931917
#> s[4]  1.00000000  0.416739030  1.4864933
#> s[5] -1.56919052 -1.462804355 -1.2383699
#> s[6]  0.54557962  0.550180667 -2.5978329
#> s[7] -0.01069215 -0.632010539 -0.1553316
#> x[1]  0.00000000  1.000000000  0.0000000
#> x[2]  0.00000000  0.000000000  1.0000000
#> x[3]  0.00000000  0.000000000  0.0000000
#> 
#> $Hsv
#> NULL
#> 
#> $nu
#> [1] 3 2 2
#> 
# check impulse response
all.equal(pseries(out$Xs, lag.max = 20), Xi)
#> [1] TRUE

Xs1 = as.stsp(Xi, method = 'echelon')
all.equal(Xs1, out$Xs)
#> [1] TRUE