1#cython: boundscheck=False 2#cython: wraparound=False 3#cython: cdivision=False 4""" 5State Space Models 6 7Author: Chad Fulton 8License: Simplified-BSD 9""" 10 11{{py: 12 13TYPES = { 14 "s": ("np.float32_t", "np.float32", "np.NPY_FLOAT32"), 15 "d": ("np.float64_t", "float", "np.NPY_FLOAT64"), 16 "c": ("np.complex64_t", "np.complex64", "np.NPY_COMPLEX64"), 17 "z": ("np.complex128_t", "complex", "np.NPY_COMPLEX128"), 18} 19 20}} 21 22# Typical imports 23import numpy as np 24import warnings 25cimport numpy as np 26cimport cython 27 28np.import_array() 29 30from statsmodels.src.math cimport * 31cimport scipy.linalg.cython_blas as blas 32cimport scipy.linalg.cython_lapack as lapack 33cimport statsmodels.tsa.statespace._tools as tools 34 35{{for prefix, types in TYPES.items()}} 36from statsmodels.tsa.statespace._initialization cimport {{prefix}}Initialization 37{{endfor}} 38 39cdef int FORTRAN = 1 40 41{{for prefix, types in TYPES.items()}} 42{{py:cython_type, dtype, typenum = types}} 43{{py: 44combined_prefix = prefix 45combined_cython_type = cython_type 46if prefix == 'c': 47 combined_prefix = 'z' 48 combined_cython_type = 'np.complex128_t' 49if prefix == 's': 50 combined_prefix = 'd' 51 combined_cython_type = 'np.float64_t' 52}} 53 54## State Space Representation 55cdef class {{prefix}}Statespace(object): 56 """ 57 {{prefix}}Statespace(obs, design, obs_intercept, obs_cov, transition, state_intercept, selection, state_cov) 58 59 *See Durbin and Koopman (2012), Chapter 4 for all notation* 60 """ 61 62 # ### State space representation 63 # 64 # $$ 65 # \begin{align} 66 # y_t & = Z_t \alpha_t + d_t + \varepsilon_t \hspace{3em} & \varepsilon_t & \sim N(0, H_t) \\\\ 67 # \alpha_{t+1} & = T_t \alpha_t + c_t + R_t \eta_t & \eta_t & \sim N(0, Q_t) \\\\ 68 # & & \alpha_1 & \sim N(a_1, P_1) 69 # \end{align} 70 # $$ 71 # 72 # $y_t$ is $p \times 1$ 73 # $\varepsilon_t$ is $p \times 1$ 74 # $\alpha_t$ is $m \times 1$ 75 # $\eta_t$ is $r \times 1$ 76 # $t = 1, \dots, T$ 77 78 # `nobs` $\equiv T$ is the length of the time-series 79 # `k_endog` $\equiv p$ is dimension of observation space 80 # `k_states` $\equiv m$ is the dimension of the state space 81 # `k_posdef` $\equiv r$ is the dimension of the state shocks 82 # *Old notation: T, n, k, g* 83 # cdef readonly int nobs, k_endog, k_states, k_posdef 84 85 # `obs` $\equiv y_t$ is the **observation vector** $(p \times T)$ 86 # `design` $\equiv Z_t$ is the **design vector** $(p \times m \times T)$ 87 # `obs_intercept` $\equiv d_t$ is the **observation intercept** $(p \times T)$ 88 # `obs_cov` $\equiv H_t$ is the **observation covariance matrix** $(p \times p \times T)$ 89 # `transition` $\equiv T_t$ is the **transition matrix** $(m \times m \times T)$ 90 # `state_intercept` $\equiv c_t$ is the **state intercept** $(m \times T)$ 91 # `selection` $\equiv R_t$ is the **selection matrix** $(m \times r \times T)$ 92 # `state_cov` $\equiv Q_t$ is the **state covariance matrix** $(r \times r \times T)$ 93 # `selected_state_cov` $\equiv R Q_t R'$ is the **selected state covariance matrix** $(m \times m \times T)$ 94 # `initial_state` $\equiv a_1$ is the **initial state mean** $(m \times 1)$ 95 # `initial_state_cov` $\equiv P_1$ is the **initial state covariance matrix** $(m \times m)$ 96 # `initial_state_cov` $\equiv P_\inf$ is the **initial diffuse state covariance matrix** $(m \times m)$ 97 # 98 # With the exception of `obs`, these are *optionally* time-varying. If they are instead time-invariant, 99 # then the dimension of length $T$ is instead of length $1$. 100 # 101 # *Note*: the initial vectors' notation 1-indexed as in Durbin and Koopman, 102 # but in the recursions below it will be 0-indexed in the Python arrays. 103 # 104 # *Old notation: y, -, mu, beta_tt_init, P_tt_init* 105 # cdef readonly {{cython_type}} [::1,:] obs, obs_intercept, state_intercept 106 # cdef readonly {{cython_type}} [:] initial_state 107 # cdef readonly {{cython_type}} [::1,:] initial_state_cov 108 # *Old notation: H, R, F, G, Q*, G Q* G'* 109 # cdef readonly {{cython_type}} [::1,:,:] design, obs_cov, transition, selection, state_cov, selected_state_cov 110 111 # `missing` is a $(p \times T)$ boolean matrix where a row is a $(p \times 1)$ vector 112 # in which the $i$th position is $1$ if $y_{i,t}$ is to be considered a missing value. 113 # *Note:* This is created as the output of np.isnan(obs). 114 # cdef readonly int [::1,:] missing 115 # `nmissing` is an `T \times 0` integer vector holding the number of *missing* observations 116 # $p - p_t$ 117 # cdef readonly int [:] nmissing 118 119 # Flag for a time-invariant model, which requires that *all* of the 120 # possibly time-varying arrays are time-invariant. 121 # cdef readonly int time_invariant 122 123 # Flag for initialization. 124 # cdef readonly int initialized 125 126 # Flags for performance improvements 127 # TODO need to add this to the UI in representation 128 # cdef public int diagonal_obs_cov 129 # cdef public int subset_design 130 # cdef public int companion_transition 131 132 # Temporary arrays 133 # cdef {{cython_type}} [::1,:] tmp 134 135 # Temporary selection arrays 136 # cdef readonly {{cython_type}} [:] selected_obs 137 # The following are contiguous memory segments which are then used to 138 # store the data in the above matrices. 139 # cdef readonly {{cython_type}} [:] selected_design 140 # cdef readonly {{cython_type}} [:] selected_obs_cov 141 142 # Temporary transformation arrays 143 # cdef readonly {{cython_type}} [::1,:] transform_cholesky 144 # cdef readonly {{cython_type}} [::1,:] transform_obs_cov 145 # cdef readonly {{cython_type}} [::1,:] transform_design 146 # cdef readonly {{cython_type}} transform_determinant 147 148 # cdef readonly {{cython_type}} [:] collapse_obs 149 # cdef readonly {{cython_type}} [:] collapse_obs_tmp 150 # cdef readonly {{cython_type}} [::1,:] collapse_design 151 # cdef readonly {{cython_type}} [::1,:] collapse_obs_cov 152 # cdef readonly {{cython_type}} [::1,:] collapse_cholesky 153 # cdef readonly {{cython_type}} collapse_loglikelihood 154 155 # Pointers 156 # cdef {{cython_type}} * _obs 157 # cdef {{cython_type}} * _design 158 # cdef {{cython_type}} * _obs_intercept 159 # cdef {{cython_type}} * _obs_cov 160 # cdef {{cython_type}} * _transition 161 # cdef {{cython_type}} * _state_intercept 162 # cdef {{cython_type}} * _selection 163 # cdef {{cython_type}} * _state_cov 164 # cdef {{cython_type}} * _selected_state_cov 165 # cdef {{cython_type}} * _initial_state 166 # cdef {{cython_type}} * _initial_state_cov 167 168 # Current location dimensions 169 # cdef int _k_endog, _k_states, _k_posdef, _k_endog2, _k_states2, _k_posdef2, _k_endogstates, _k_statesposdef 170 # cdef int _nmissing 171 172 # ### Initialize state space model 173 # *Note*: The initial state and state covariance matrix must be provided. 174 def __init__(self, 175 {{cython_type}} [::1,:] obs, 176 {{cython_type}} [::1,:,:] design, 177 {{cython_type}} [::1,:] obs_intercept, 178 {{cython_type}} [::1,:,:] obs_cov, 179 {{cython_type}} [::1,:,:] transition, 180 {{cython_type}} [::1,:] state_intercept, 181 {{cython_type}} [::1,:,:] selection, 182 {{cython_type}} [::1,:,:] state_cov, 183 diagonal_obs_cov=-1): 184 185 # Local variables 186 cdef: 187 int t, i, j 188 np.npy_intp dim1[1] 189 np.npy_intp dim2[2] 190 np.npy_intp dim3[3] 191 192 # #### State space representation variables 193 # **Note**: these arrays share data with the versions defined in 194 # Python and passed to this constructor, so if they are updated in 195 # Python they will also be updated here. 196 self.obs = obs 197 self.design = design 198 self.obs_intercept = obs_intercept 199 self.obs_cov = obs_cov 200 self.transition = transition 201 self.state_intercept = state_intercept 202 self.selection = selection 203 self.state_cov = state_cov 204 205 # Dimensions 206 self.k_endog = obs.shape[0] 207 self.k_states = selection.shape[0] 208 self.k_posdef = selection.shape[1] 209 self.nobs = obs.shape[1] 210 211 # #### Validate matrix dimensions 212 # 213 # Make sure that the given state-space matrices have consistent sizes 214 tools.validate_matrix_shape('design', &self.design.shape[0], 215 self.k_endog, self.k_states, self.nobs) 216 tools.validate_vector_shape('observation intercept', &self.obs_intercept.shape[0], 217 self.k_endog, self.nobs) 218 tools.validate_matrix_shape('observation covariance matrix', &self.obs_cov.shape[0], 219 self.k_endog, self.k_endog, self.nobs) 220 tools.validate_matrix_shape('transition', &self.transition.shape[0], 221 self.k_states, self.k_states, self.nobs) 222 tools.validate_vector_shape('state intercept', &self.state_intercept.shape[0], 223 self.k_states, self.nobs) 224 tools.validate_matrix_shape('state covariance matrix', &self.state_cov.shape[0], 225 self.k_posdef, self.k_posdef, self.nobs) 226 227 # Check for a time-invariant model 228 self.time_invariant = ( 229 self.design.shape[2] == 1 and 230 self.obs_intercept.shape[1] == 1 and 231 self.obs_cov.shape[2] == 1 and 232 self.transition.shape[2] == 1 and 233 self.state_intercept.shape[1] == 1 and 234 self.selection.shape[2] == 1 and 235 self.state_cov.shape[2] == 1 236 ) 237 238 # Set the flags for initialization to be false 239 self.initialized = False 240 self.initialized_diffuse = False 241 self.initialized_stationary = False 242 243 self.diagonal_obs_cov = diagonal_obs_cov 244 self._diagonal_obs_cov = -1 245 246 # Allocate selected state covariance matrix 247 dim3[0] = self.k_states; dim3[1] = self.k_states; dim3[2] = 1; 248 # (we only allocate memory for time-varying array if necessary) 249 if self.state_cov.shape[2] > 1 or self.selection.shape[2] > 1: 250 dim3[2] = self.nobs 251 self.selected_state_cov = np.PyArray_ZEROS(3, dim3, {{typenum}}, FORTRAN) 252 253 # Handle missing data 254 self.missing = np.array(np.isnan(obs), dtype=np.int32, order="F") 255 self.nmissing = np.array(np.sum(self.missing, axis=0), dtype=np.int32) 256 self.has_missing = np.sum(self.nmissing) > 0 257 258 # Create the temporary array 259 # Holds arrays of dimension $(m \times m)$ 260 dim2[0] = self.k_states; dim2[1] = max(self.k_states, self.k_posdef); 261 self.tmp = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN) 262 263 # Arrays for initialization 264 dim1[0] = self.k_states; 265 self.initial_state = np.PyArray_ZEROS(1, dim1, {{typenum}}, FORTRAN) 266 dim2[0] = self.k_states; dim2[1] = self.k_states; 267 self.initial_state_cov = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN) 268 dim2[0] = self.k_states; dim2[1] = self.k_states; 269 self.initial_diffuse_state_cov = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN) 270 271 # Arrays for missing data 272 dim1[0] = self.k_endog; 273 self.selected_obs = np.PyArray_ZEROS(1, dim1, {{typenum}}, FORTRAN) 274 dim1[0] = self.k_endog; 275 self.selected_obs_intercept = np.PyArray_ZEROS(1, dim1, {{typenum}}, FORTRAN) 276 dim1[0] = self.k_endog * self.k_states; 277 self.selected_design = np.PyArray_ZEROS(1, dim1, {{typenum}}, FORTRAN) 278 dim1[0] = self.k_endog**2; 279 self.selected_obs_cov = np.PyArray_ZEROS(1, dim1, {{typenum}}, FORTRAN) 280 281 # Arrays for transformations 282 dim2[0] = self.k_endog; dim2[1] = self.k_endog; 283 self.transform_cholesky = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN) 284 self.transform_obs_cov = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN) 285 dim2[0] = self.k_endog; dim2[1] = self.k_states; 286 self.transform_design = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN) 287 dim1[0] = self.k_endog; 288 self.transform_obs_intercept = np.PyArray_ZEROS(1, dim1, {{typenum}}, FORTRAN) 289 290 dim1[0] = self.k_states; 291 self.collapse_obs = np.PyArray_ZEROS(1, dim1, {{typenum}}, FORTRAN) 292 self.collapse_obs_tmp = np.PyArray_ZEROS(1, dim1, {{typenum}}, FORTRAN) 293 dim2[0] = self.k_states; dim2[1] = self.k_states; 294 self.collapse_design = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN) 295 self.collapse_obs_cov = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN) 296 self.collapse_cholesky = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN) 297 298 # Initialize location 299 self.t = 0 300 self._previous_t = 0 301 302 # Initialize dimensions 303 self.set_dimensions(self.k_endog, self.k_states, self.k_posdef) 304 305 def __reduce__(self): 306 init = (np.array(self.obs, copy=True, order='F'), np.array(self.design, copy=True, order='F'), 307 np.array(self.obs_intercept, copy=True, order='F'), np.array(self.obs_cov, copy=True, order='F'), 308 np.array(self.transition, copy=True, order='F'), np.array(self.state_intercept, copy=True, order='F'), 309 np.array(self.selection, copy=True, order='F'), np.array(self.state_cov, copy=True, order='F'), 310 self.diagonal_obs_cov) 311 state = {'initialized': self.initialized, 312 'initialized_diffuse': self.initialized_diffuse, 313 'initialized_stationary': self.initialized_stationary, 314 'initial_state': None, 315 'initial_state_cov': None, 316 'initial_diffuse_state_cov': None, 317 'missing': np.array(self.missing, copy=True, order='F'), 318 'nmissing': np.array(self.nmissing, copy=True, order='F'), 319 'has_missing': self.has_missing, 320 'tmp': np.array(self.tmp, copy=True, order='F'), 321 'selected_state_cov': np.array(self.selected_state_cov, copy=True, order='F'), 322 'selected_obs': np.array(self.selected_obs, copy=True, order='F'), 323 'selected_obs_intercept': np.array(self.selected_obs_intercept, copy=True, order='F'), 324 'selected_design': np.array(self.selected_design, copy=True, order='F'), 325 'selected_obs_cov': np.array(self.selected_obs_cov, copy=True, order='F'), 326 'transform_cholesky': np.array(self.transform_cholesky, copy=True, order='F'), 327 'transform_obs_cov': np.array(self.transform_obs_cov, copy=True, order='F'), 328 'transform_design': np.array(self.transform_design, copy=True, order='F'), 329 'collapse_obs': np.array(self.collapse_obs, copy=True, order='F'), 330 'collapse_obs_tmp': np.array(self.collapse_obs_tmp, copy=True, order='F'), 331 'collapse_design': np.array(self.collapse_design, copy=True, order='F'), 332 'collapse_obs_cov': np.array(self.collapse_obs_cov, copy=True, order='F'), 333 'collapse_cholesky': np.array(self.collapse_cholesky, copy=True, order='F'), 334 't': self.t, 335 'collapse_loglikelihood': self.collapse_loglikelihood, 336 'companion_transition': self.companion_transition, 337 'transform_determinant': self.transform_determinant, 338 } 339 if self.initialized: 340 state['initial_state'] = np.array(self.initial_state, copy=True, order='F') 341 state['initial_state_cov'] = np.array(self.initial_state_cov, copy=True, order='F') 342 state['initial_diffuse_state_cov'] = np.array(self.initial_diffuse_state_cov, copy=True, order='F') 343 return (self.__class__, init, state) 344 345 def __setstate__(self, state): 346 self.initial_state = state['initial_state'] 347 self.initial_state_cov = state['initial_state_cov'] 348 self.initial_diffuse_state_cov = state['initial_diffuse_state_cov'] 349 self.initialized = state['initialized'] 350 self.initialized_diffuse = state['initialized_diffuse'] 351 self.initialized_stationary = state['initialized_stationary'] 352 self.selected_state_cov = state['selected_state_cov'] 353 self.missing = state['missing'] 354 self.nmissing =state['nmissing'] 355 self.has_missing = state['has_missing'] 356 self.tmp = state['tmp'] 357 self.selected_obs = state['selected_obs'] 358 self.selected_obs_intercept = state['selected_obs_intercept'] 359 self.selected_design = state['selected_design'] 360 self.selected_obs_cov =state['selected_obs_cov'] 361 self.transform_cholesky = state['transform_cholesky'] 362 self.transform_obs_cov = state['transform_obs_cov'] 363 self.transform_design = state['transform_design'] 364 self.collapse_obs = state['collapse_obs'] 365 self.collapse_obs_tmp = state['collapse_obs_tmp'] 366 self.collapse_design = state['collapse_design'] 367 self.collapse_obs_cov = state['collapse_obs_cov'] 368 self.collapse_cholesky = state['collapse_cholesky'] 369 self.t = state['t'] 370 self.collapse_loglikelihood = state['collapse_loglikelihood'] 371 self.companion_transition = state['companion_transition'] 372 self.transform_determinant = state['transform_determinant'] 373 374 def initialize(self, init, offset=0, complex_step=False, clear=True): 375 cdef {{prefix}}Initialization _init 376 # Clear initial arrays 377 if clear: 378 self.initial_state[:] = 0 379 self.initial_diffuse_state_cov[:] = 0 380 self.initial_state_cov[:] = 0 381 382 # If using global initialization, compute the actual elements and 383 # return them 384 self.initialized_diffuse = False 385 self.initialized_stationary = False 386 if init.initialization_type is not None: 387 init._initialize_initialization(prefix='{{prefix}}') 388 _init = init._initializations['{{prefix}}'] 389 _init.initialize(init.initialization_type, offset, self, 390 self.initial_state, 391 self.initial_diffuse_state_cov, 392 self.initial_state_cov, complex_step) 393 if init.initialization_type == 'diffuse': 394 self.initialized_diffuse = True 395 if init.initialization_type == 'stationary': 396 self.initialized_stationary = True 397 # Otherwise, if using blocks, initialize each of the blocks 398 else: 399 for block_index, block_init in init.blocks.items(): 400 self.initialize(block_init, offset=offset + block_index[0], 401 complex_step=complex_step, clear=False) 402 403 if not self.initialized: 404 self.initialized = True 405 406 # ## Initialize: known values 407 # 408 # Initialize the filter with specific values, assumed to be known with 409 # certainty or else as filled with parameters from a maximum likelihood 410 # estimation run. 411 def initialize_known(self, {{cython_type}} [:] initial_state, {{cython_type}} [::1,:] initial_state_cov): 412 """ 413 initialize_known(initial_state, initial_state_cov) 414 """ 415 tools.validate_vector_shape('initial state', &initial_state.shape[0], self.k_states, None) 416 tools.validate_matrix_shape('initial state covariance', &initial_state_cov.shape[0], self.k_states, self.k_states, None) 417 418 self.initial_state = initial_state 419 self.initial_state_cov = initial_state_cov 420 self.initial_diffuse_state_cov[:] = 0 421 422 self.initialized = True 423 424 # ## Initialize: approximate diffuse priors 425 # 426 # Durbin and Koopman note that this initialization should only be coupled 427 # with the standard Kalman filter for "approximate exploratory work" and 428 # can lead to "large rounding errors" (p. 125). 429 # 430 # *Note:* see Durbin and Koopman section 5.6.1 431 def initialize_approximate_diffuse(self, {{cython_type}} variance=1e2): 432 """ 433 initialize_approximate_diffuse(variance=1e2) 434 """ 435 cdef np.npy_intp dim[1] 436 dim[0] = self.k_states 437 self.initial_state = np.PyArray_ZEROS(1, dim, {{typenum}}, FORTRAN) 438 self.initial_state_cov = np.eye(self.k_states, dtype={{dtype}}).T * variance 439 self.initial_diffuse_state_cov[:] = 0 440 441 self.initialized = True 442 443 # ## Initialize: stationary process 444 # *Note:* see Durbin and Koopman section 5.6.2 445 def initialize_stationary(self, complex_step=False): 446 """ 447 initialize_stationary() 448 """ 449 cdef np.npy_intp dim1[1] 450 cdef np.npy_intp dim2[2] 451 cdef int i, info, inc = 1 452 cdef int k_states2 = self.k_states**2 453 cdef np.float64_t asum, tol = 1e-9 454 cdef {{cython_type}} scalar 455 cdef int [::1,:] ipiv 456 457 # Create selected state covariance matrix 458 {{prefix}}select_cov(self.k_states, self.k_posdef, 459 &self.tmp[0,0], 460 &self.selection[0,0,0], 461 &self.state_cov[0,0,0], 462 &self.selected_state_cov[0,0,0]) 463 464 # Initial state mean 465 {{if combined_prefix == 'd'}} 466 asum = blas.{{prefix}}asum(&self.k_states, &self.state_intercept[0, 0], &inc) 467 {{elif prefix == 'c'}} 468 asum = blas.scasum(&self.k_states, &self.state_intercept[0, 0], &inc) 469 {{else}} 470 asum = blas.dzasum(&self.k_states, &self.state_intercept[0, 0], &inc) 471 {{endif}} 472 473 dim1[0] = self.k_states 474 self.initial_state = np.PyArray_ZEROS(1, dim1, {{typenum}}, FORTRAN) 475 if asum > tol: 476 dim2[0] = self.k_states 477 dim2[1] = self.k_states 478 ipiv = np.PyArray_ZEROS(2, dim2, np.NPY_INT32, FORTRAN) 479 480 # I - T 481 blas.{{prefix}}copy(&k_states2, &self.transition[0,0,0], &inc, 482 &self.tmp[0,0], &inc) 483 scalar = -1.0 484 blas.{{prefix}}scal(&k_states2, &scalar, &self.tmp[0, 0], &inc) 485 for i in range(self.k_states): 486 self.tmp[i, i] = self.tmp[i, i] + 1 487 488 # c 489 blas.{{prefix}}copy(&self.k_states, &self.state_intercept[0,0], &inc, 490 &self.initial_state[0], &inc) 491 492 # Solve (I - T) x = c 493 lapack.{{prefix}}getrf(&self.k_states, &self.k_states, &self.tmp[0, 0], &self.k_states, 494 &ipiv[0, 0], &info) 495 lapack.{{prefix}}getrs('N', &self.k_states, &inc, &self.tmp[0, 0], &self.k_states, 496 &ipiv[0, 0], &self.initial_state[0], &self.k_states, &info) 497 498 dim2[0] = self.k_states; dim2[1] = self.k_states; 499 self.initial_state_cov = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN) 500 501 # Create a copy of the transition matrix (to avoid overwriting it) 502 blas.{{prefix}}copy(&k_states2, &self.transition[0,0,0], &inc, 503 &self.tmp[0,0], &inc) 504 505 # Copy the selected state covariance to the initial state covariance 506 # (it will be overwritten with the appropriate matrix) 507 blas.{{prefix}}copy(&k_states2, &self.selected_state_cov[0,0,0], &inc, 508 &self.initial_state_cov[0,0], &inc) 509 510 # Solve the discrete Lyapunov equation to the get initial state 511 # covariance matrix 512 tools._{{prefix}}solve_discrete_lyapunov(&self.tmp[0,0], &self.initial_state_cov[0,0], self.k_states, complex_step) 513 514 self.initial_diffuse_state_cov[:] = 0 515 516 self.initialized = True 517 518 def __iter__(self): 519 return self 520 521 def __next__(self): 522 """ 523 Advance to the next location 524 """ 525 if self.t >= self.nobs: 526 raise StopIteration 527 else: 528 self.seek(self.t+1, 0, 0) 529 530 cpdef seek(self, unsigned int t, unsigned int transform_diagonalize, unsigned int transform_generalized_collapse, unsigned int reset=False): 531 self._previous_t = self.t 532 533 # Set the global time indicator, if valid 534 if t >= self.nobs: 535 raise IndexError("Observation index out of range") 536 self.t = t 537 538 # Indices for possibly time-varying arrays 539 cdef: 540 int k_endog 541 int design_t = 0 542 int obs_intercept_t = 0 543 int obs_cov_t = 0 544 int transition_t = 0 545 int state_intercept_t = 0 546 int selection_t = 0 547 int state_cov_t = 0 548 549 # Get indices for possibly time-varying arrays 550 if not self.time_invariant: 551 if self.design.shape[2] > 1: design_t = t 552 if self.obs_intercept.shape[1] > 1: obs_intercept_t = t 553 if self.obs_cov.shape[2] > 1: obs_cov_t = t 554 if self.transition.shape[2] > 1: transition_t = t 555 if self.state_intercept.shape[1] > 1: state_intercept_t = t 556 if self.selection.shape[2] > 1: selection_t = t 557 if self.state_cov.shape[2] > 1: state_cov_t = t 558 559 # Initialize object-level pointers to statespace arrays 560 self._obs = &self.obs[0, t] 561 self._design = &self.design[0, 0, design_t] 562 self._obs_intercept = &self.obs_intercept[0, obs_intercept_t] 563 self._obs_cov = &self.obs_cov[0, 0, obs_cov_t] 564 self._transition = &self.transition[0, 0, transition_t] 565 self._state_intercept = &self.state_intercept[0, state_intercept_t] 566 self._selection = &self.selection[0, 0, selection_t] 567 self._state_cov = &self.state_cov[0, 0, state_cov_t] 568 569 # Initialize object-level pointers to initialization 570 if not self.initialized: 571 raise RuntimeError("Statespace model not initialized.") 572 self._initial_state = &self.initial_state[0] 573 self._initial_state_cov = &self.initial_state_cov[0,0] 574 self._initial_diffuse_state_cov = &self.initial_diffuse_state_cov[0,0] 575 576 # Create the selected state covariance matrix 577 self.select_state_cov(t) 578 579 # Handle missing data 580 # Note: this modifies object pointers and _* dimensions 581 k_endog = self.select_missing(t) 582 583 # Set dimensions 584 self.set_dimensions(k_endog, self.k_states, self.k_posdef) 585 586 # Handle transformations 587 self.transform(t, self._previous_t, transform_diagonalize, transform_generalized_collapse, reset) 588 589 cdef void set_dimensions(self, unsigned int k_endog, unsigned int k_states, unsigned int k_posdef): 590 self._k_endog = k_endog 591 self._k_states = k_states 592 self._k_posdef = k_posdef 593 self._k_endog2 = k_endog**2 594 self._k_states2 = k_states**2 595 self._k_posdef2 = k_posdef**2 596 self._k_endogstates = k_endog * k_states 597 self._k_statesposdef = k_states * k_posdef 598 599 cdef void select_state_cov(self, unsigned int t): 600 cdef int selected_state_cov_t = 0 601 602 # ### Get selected state covariance matrix 603 if t == 0 or self.selected_state_cov.shape[2] > 1: 604 selected_state_cov_t = t 605 self._selected_state_cov = &self.selected_state_cov[0, 0, selected_state_cov_t] 606 607 {{prefix}}select_cov(self.k_states, self.k_posdef, 608 &self.tmp[0,0], 609 self._selection, 610 self._state_cov, 611 self._selected_state_cov) 612 else: 613 self._selected_state_cov = &self.selected_state_cov[0, 0, 0] 614 615 cdef int select_missing(self, unsigned int t): 616 # Note: this assumes that object pointers are already initialized 617 # Note: this assumes that transform_... will be done *later* 618 cdef int k_endog = self.k_endog 619 620 # Set the current iteration nmissing 621 self._nmissing = self.nmissing[t] 622 623 # ### Perform missing selections 624 # In Durbin and Koopman (2012), these are represented as matrix 625 # multiplications, i.e. $Z_t^* = W_t Z_t$ where $W_t$ is a row 626 # selection matrix (it contains a subset of rows of the identity 627 # matrix). 628 # 629 # It's more efficient, though, to just copy over the data directly, 630 # which is what is done here. Note that the `selected_*` arrays are 631 # defined as single-dimensional, so the assignment indexes below are 632 # set such that the arrays can be interpreted by the BLAS and LAPACK 633 # functions as two-dimensional, column-major arrays. 634 # 635 # In the case that all data is missing (e.g. this is what happens in 636 # forecasting), we actually set do not change the dimension, but we set 637 # the design matrix to the zeros array. 638 if self._nmissing == self.k_endog: 639 self._select_missing_entire_obs(t) 640 elif self._nmissing > 0: 641 self._select_missing_partial_obs(t) 642 k_endog = self.k_endog - self._nmissing 643 644 # Return the number of non-missing endogenous variables 645 return k_endog 646 647 cdef void _select_missing_entire_obs(self, unsigned int t): 648 cdef: 649 int i, j 650 651 # Design matrix is set to zeros 652 for i in range(self.k_states): 653 for j in range(self.k_endog): 654 self.selected_design[j + i*self.k_endog] = 0.0 655 self._design = &self.selected_design[0] 656 657 cdef void _select_missing_partial_obs(self, unsigned int t): 658 cdef: 659 int i, j, k, l 660 int inc = 1 661 int design_t = 0 662 int obs_cov_t = 0 663 int k_endog = self.k_endog - self._nmissing 664 665 k = 0 666 for i in range(self.k_endog): 667 if not self.missing[i, t]: 668 669 self.selected_obs[k] = self._obs[i] 670 self.selected_obs_intercept[k] = self._obs_intercept[i] 671 672 # i is rows, k is rows 673 blas.{{prefix}}copy(&self.k_states, 674 &self._design[i], &self.k_endog, 675 &self.selected_design[k], &k_endog) 676 677 # i, k is columns, j, l is rows 678 l = 0 679 for j in range(self.k_endog): 680 if not self.missing[j, t]: 681 self.selected_obs_cov[l + k*k_endog] = self._obs_cov[j + i*self.k_endog] 682 l += 1 683 k += 1 684 self._obs = &self.selected_obs[0] 685 self._obs_intercept = &self.selected_obs_intercept[0] 686 self._design = &self.selected_design[0] 687 self._obs_cov = &self.selected_obs_cov[0] 688 689 cdef void transform(self, unsigned int t, unsigned int previous_t, unsigned int transform_diagonalize, unsigned int transform_generalized_collapse, unsigned int reset=False) except *: 690 # Reset the collapsed loglikelihood 691 self.collapse_loglikelihood = 0 692 693 if transform_generalized_collapse and not self._k_endog <= self._k_states: 694 k_endog = self.transform_generalized_collapse(t, previous_t, reset) 695 # Reset dimensions 696 self.set_dimensions(k_endog, self._k_states, self._k_posdef) 697 elif transform_diagonalize and not (self.diagonal_obs_cov == 1): 698 self.transform_diagonalize(t, previous_t, reset) 699 700 cdef void transform_diagonalize(self, unsigned int t, unsigned int previous_t, unsigned int reset=False) except *: 701 # Note: this assumes that initialize_object_pointers has *already* been done 702 # Note: this assumes that select_missing has *already* been done 703 # TODO need unit tests, especially for the missing case 704 # TODO need to also transform observation intercept 705 cdef: 706 int i, j, inc=1 707 int obs_cov_t = 0 708 int info 709 int reset_missing 710 int diagonal_obs_cov 711 {{cython_type}} * _transform_obs_cov = &self.transform_obs_cov[0, 0] 712 {{cython_type}} * _transform_cholesky = &self.transform_cholesky[0, 0] 713 714 if self.obs_cov.shape[2] > 1: 715 obs_cov_t = t 716 717 # Handle missing data 718 if self.nmissing[t] == self.k_endog: 719 return 720 reset_missing = 0 721 for i in range(self.k_endog): 722 reset_missing = reset_missing + (not self.missing[i,t] == self.missing[i,previous_t]) 723 724 # If the flag for a diagonal covariancem matrix is not set globally in 725 # the model one way or the other, then we need to check 726 diagonal_obs_cov = self.diagonal_obs_cov 727 if diagonal_obs_cov == -1: 728 # We don't need to check for a diagonal covariance matrix each t, 729 # except in the following cases: 730 if self._diagonal_obs_cov == -1 or t == 0 or self.obs_cov.shape[2] > 1 or reset_missing or reset: 731 diagonal_obs_cov = 1 732 for i in range(self.k_endog): 733 for j in range(self.k_endog): 734 if i == j: 735 continue 736 if not ({{combined_prefix}}abs(self.obs_cov[i, j, obs_cov_t]) < 1e-9): 737 diagonal_obs_cov = 0 738 break 739 # Otherwise, we use whatever value was produced last period 740 else: 741 diagonal_obs_cov = self._diagonal_obs_cov 742 self._diagonal_obs_cov = diagonal_obs_cov 743 if diagonal_obs_cov == 1: 744 return 745 746 # If we have a non-diagonal obs cov, we need to compute the cholesky 747 # decomposition of *self._obs_cov 748 if t == 0 or self.obs_cov.shape[2] > 1 or reset_missing or reset: 749 # LDL decomposition 750 blas.{{prefix}}copy(&self._k_endog2, self._obs_cov, &inc, _transform_cholesky, &inc) 751 info = tools._{{prefix}}ldl(_transform_cholesky, self._k_endog) 752 753 # Check for errors 754 if info > 0: 755 warnings.warn("Positive semi-definite observation covariance matrix encountered at period %d" % t) 756 elif info < 0: 757 raise np.linalg.LinAlgError('Invalid value in LDL factorization of observation covariance matrix encountered at period %d' % t) 758 759 # Currently both L and D are stored in transform_cholesky 760 for i in range(self._k_endog): # i is rows 761 for j in range(self._k_endog): # j is columns 762 # Diagonal elements come from the diagonal 763 if i == j: 764 _transform_obs_cov[i + i * self._k_endog] = _transform_cholesky[i + i * self._k_endog] 765 # Other elements are zero 766 else: 767 _transform_obs_cov[i + j * self._k_endog] = 0 768 769 # Zero out the upper triangle of the cholesky 770 if j > i: 771 _transform_cholesky[i + j * self._k_endog] = 0 772 773 # Convert from L to C simply by setting the diagonal elements to ones 774 _transform_cholesky[i + i * self._k_endog] = 1 775 776 # Solve for y_t^* 777 # (unless this is a completely missing observation) 778 # TODO: note that this can cause problems if this function is run twice 779 if not self._nmissing == self.k_endog: 780 # If we have some missing elements, selected_obs is already populated 781 if self._nmissing == 0: 782 blas.{{prefix}}copy(&self._k_endog, &self.obs[0,t], &inc, &self.selected_obs[0], &inc) 783 lapack.{{prefix}}trtrs("L", "N", "U", &self._k_endog, &inc, 784 _transform_cholesky, &self.k_endog, 785 &self.selected_obs[0], &self._k_endog, &info) 786 787 # Check for errors 788 if info > 0: 789 raise np.linalg.LinAlgError('Singular factorization of observation covariance matrix encountered at period %d' % t) 790 elif info < 0: 791 raise np.linalg.LinAlgError('Invalid value in factorization of observation covariance matrix encountered at period %d' % t) 792 793 # Setup the pointer 794 self._obs = &self.selected_obs[0] 795 796 # Solve for d_t^*, if necessary 797 if t == 0 or self.obs_intercept.shape[1] > 1 or reset_missing or reset: 798 blas.{{prefix}}copy(&self._k_endog, self._obs_intercept, &inc, &self.transform_obs_intercept[0], &inc) 799 lapack.{{prefix}}trtrs("L", "N", "U", &self._k_endog, &inc, 800 _transform_cholesky, &self.k_endog, 801 &self.transform_obs_intercept[0], &self._k_endog, 802 &info) 803 804 # Solve for Z_t^*, if necessary 805 if t == 0 or self.design.shape[2] > 1 or reset_missing or reset: 806 blas.{{prefix}}copy(&self._k_endogstates, self._design, &inc, &self.transform_design[0,0], &inc) 807 lapack.{{prefix}}trtrs("L", "N", "U", &self._k_endog, &self._k_states, 808 _transform_cholesky, &self.k_endog, 809 &self.transform_design[0,0], &self._k_endog, 810 &info) 811 812 # Check for errors 813 if info > 0: 814 raise np.linalg.LinAlgError('Singular factorization of observation covariance matrix encountered at period %d' % t) 815 elif info < 0: 816 raise np.linalg.LinAlgError('Invalid value in factorization of observation covariance matrix encountered at period %d' % t) 817 818 # Setup final pointers 819 self._design = &self.transform_design[0,0] 820 self._obs_cov = &self.transform_obs_cov[0,0] 821 self._obs_intercept = &self.transform_obs_intercept[0] 822 823 cdef int transform_generalized_collapse(self, unsigned int t, unsigned int previous_t, unsigned int reset=True) except *: 824 # Note: this assumes that initialize_object_pointers has *already* been done 825 # Note: this assumes that select_missing has *already* been done 826 # TODO need unit tests, especially for the missing case 827 cdef: 828 int i, j, k, l, inc=1 829 int obs_cov_t, design_t 830 int info 831 int reset_missing 832 {{cython_type}} alpha = 1.0 833 {{cython_type}} beta = 0.0 834 {{cython_type}} gamma = -1.0 835 int k_states = self._k_states 836 int k_states2 = self._k_states2 837 int k_endogstates = self._k_endogstates 838 839 # $y_t^* = \bar A^* y_t = C_t Z_t' H_t^{-1} y_t$ 840 # $Z_t^* = C_t^{-1}$ 841 # $H_t^* = I_m$ 842 843 # Make sure we have enough observations to perform collapse 844 if self.k_endog < self.k_states: 845 raise RuntimeError('Cannot collapse observation vector it the' 846 ' state dimension is larger than the dimension' 847 ' of the observation vector.') 848 849 # Adjust for a VAR transition (i.e. design = [#, 0], where the zeros 850 # correspond to all states except the first k_posdef states) 851 if self.subset_design: 852 k_states = self._k_posdef 853 k_states2 = self._k_posdef2 854 k_endogstates = self._k_endog * self._k_posdef 855 856 # Handle missing data 857 if self.nmissing[t] == self.k_endog: 858 return self.k_states 859 reset_missing = 0 860 for i in range(self.k_endog): 861 reset_missing = reset_missing + (not self.missing[i,t] == self.missing[i,previous_t]) 862 863 # Initialize the transformation 864 if self.collapse_obs_cov[0,0] == 0: 865 # Set H_t^* to identity 866 for i in range(k_states): 867 self.collapse_obs_cov[i,i] = 1 868 869 # Make sure we do not have an observation intercept 870 if not np.sum(self.obs_intercept) == 0 or self.obs_intercept.shape[2] > 1: 871 raise RuntimeError('The observation collapse transformation' 872 ' does not currently support an observation' 873 ' intercept.') 874 875 # Perform the Cholesky decomposition of H_t, if necessary 876 if t == 0 or self.obs_cov.shape[2] > 1 or reset_missing or reset: 877 # Cholesky decomposition: $H = L L'$ 878 blas.{{prefix}}copy(&self._k_endog2, self._obs_cov, &inc, &self.transform_cholesky[0,0], &inc) 879 lapack.{{prefix}}potrf("L", &self._k_endog, &self.transform_cholesky[0,0], &self._k_endog, &info) 880 881 # Check for errors 882 if info > 0: 883 raise np.linalg.LinAlgError('Non-positive-definite observation covariance matrix encountered at period %d' % t) 884 elif info < 0: 885 raise np.linalg.LinAlgError('Invalid value in observation covariance matrix encountered at period %d' % t) 886 887 # Calculate the determinant (just the squared product of the 888 # diagonals, in the Cholesky decomposition case) 889 self.transform_determinant = 0. 890 for i in range(self._k_endog): 891 j = i * (self._k_endog + 1) 892 k = j % self.k_endog 893 l = j // self.k_endog 894 if not self.transform_cholesky[k, l] == 0: 895 self.transform_determinant = self.transform_determinant + {{combined_prefix}}log(self.transform_cholesky[k, l]) 896 self.transform_determinant = 2 * self.transform_determinant 897 898 # Get $Z_t \equiv C^{-1}$, if necessary 899 if t == 0 or self.obs_cov.shape[2] > 1 or self.design.shape[2] > 1 or reset_missing or reset: 900 # Calculate $H_t^{-1} Z_t \equiv (Z_t' H_t^{-1})'$ via Cholesky solver 901 blas.{{prefix}}copy(&self._k_endogstates, self._design, &inc, &self.transform_design[0,0], &inc) 902 lapack.{{prefix}}potrs("L", &self._k_endog, &k_states, 903 &self.transform_cholesky[0,0], &self._k_endog, 904 &self.transform_design[0,0], &self._k_endog, 905 &info) 906 907 # Check for errors 908 if not info == 0: 909 raise np.linalg.LinAlgError('Invalid value in calculation of H_t^{-1}Z matrix encountered at period %d' % t) 910 911 # Calculate $(H_t^{-1} Z_t)' Z_t$ 912 # $(m \times m) = (m \times p) (p \times p) (p \times m)$ 913 blas.{{prefix}}gemm("T", "N", &k_states, &k_states, &self._k_endog, 914 &alpha, self._design, &self._k_endog, 915 &self.transform_design[0,0], &self._k_endog, 916 &beta, &self.collapse_cholesky[0,0], &self._k_states) 917 918 # Calculate $(Z_t' H_t^{-1} Z_t)^{-1}$ via Cholesky inversion 919 lapack.{{prefix}}potrf("U", &k_states, &self.collapse_cholesky[0,0], &self.k_states, &info) 920 # Check for errors 921 if info > 0: 922 raise np.linalg.LinAlgError('Non-positive-definite ZHZ matrix encountered at period %d' % t) 923 elif info < 0: 924 raise np.linalg.LinAlgError('Invalid value in ZHZ matrix encountered at period %d' % t) 925 lapack.{{prefix}}potri("U", &k_states, &self.collapse_cholesky[0,0], &self.k_states, &info) 926 927 # Calculate $C_t$ (the upper triangular cholesky decomposition of $(Z_t' H_t^{-1} Z_t)^{-1}$) 928 lapack.{{prefix}}potrf("U", &k_states, &self.collapse_cholesky[0,0], &self.k_states, &info) 929 930 # Check for errors 931 if info > 0: 932 raise np.linalg.LinAlgError('Non-positive-definite C matrix encountered at period %d' % t) 933 elif info < 0: 934 raise np.linalg.LinAlgError('Invalid value in C matrix encountered at period %d' % t) 935 936 # Calculate $C_t'^{-1} \equiv Z_t$ 937 # Do so by solving the system: $C_t' x = I$ 938 # (Recall that collapse_obs_cov is an identity matrix) 939 blas.{{prefix}}copy(&self._k_states2, &self.collapse_obs_cov[0,0], &inc, &self.collapse_design[0,0], &inc) 940 lapack.{{prefix}}trtrs("U", "T", "N", &k_states, &k_states, 941 &self.collapse_cholesky[0,0], &self._k_states, 942 &self.collapse_design[0,0], &self._k_states, 943 &info) 944 945 # Calculate $\bar y_t^* = \bar A_t^* y_t = C_t Z_t' H_t^{-1} y_t$ 946 # (unless this is a completely missing observation) 947 self.collapse_loglikelihood = 0 948 if not self._nmissing == self.k_endog: 949 # If we have some missing elements, selected_obs is already populated 950 if self._nmissing == 0: 951 blas.{{prefix}}copy(&self.k_endog, &self.obs[0,t], &inc, &self.selected_obs[0], &inc) 952 # $\\# = Z_t' H_t^{-1} y_t$ 953 blas.{{prefix}}gemv("T", &self._k_endog, &k_states, 954 &alpha, &self.transform_design[0,0], &self._k_endog, 955 &self.selected_obs[0], &inc, 956 &beta, &self.collapse_obs[0], &inc) 957 # $y_t^* = C_t \\#$ 958 blas.{{prefix}}trmv("U", "N", "N", &k_states, 959 &self.collapse_cholesky[0,0], &self._k_states, 960 &self.collapse_obs[0], &inc) 961 962 # Get residuals for loglikelihood calculation 963 # Note: Durbin and Koopman (2012) appears to have an error in the 964 # formula here. They have $e_t = y_t - Z_t \bar y_t^*$, whereas it 965 # should be: $e_t = y_t - Z_t C_t' \bar y_t^*$ 966 # See Jungbacker and Koopman (2014), section 2.5 where $e_t$ is 967 # defined. In this case, $Z_t^dagger = Z_t C_t$ where 968 # $C_t C_t' = (Z_t' \Sigma_\varepsilon^{-1} Z_t)^{-1}$. 969 # 970 971 # $ \\# = C_t' y_t^*$ 972 blas.{{prefix}}copy(&k_states, &self.collapse_obs[0], &inc, &self.collapse_obs_tmp[0], &inc) 973 blas.{{prefix}}trmv("U", "T", "N", &k_states, 974 &self.collapse_cholesky[0,0], &self._k_states, 975 &self.collapse_obs_tmp[0], &inc) 976 977 # $e_t = - Z_t C_t' y_t^* + y_t$ 978 blas.{{prefix}}gemv("N", &self._k_endog, &k_states, 979 &gamma, self._design, &self._k_endog, 980 &self.collapse_obs_tmp[0], &inc, 981 &alpha, &self.selected_obs[0], &inc) 982 983 # Calculate e_t' H_t^{-1} e_t via Cholesky solver 984 # $H_t^{-1} = (L L')^{-1} = L^{-1}' L^{-1}$ 985 # So we want $e_t' L^{-1}' L^{-1} e_t = (L^{-1} e_t)' L^{-1} e_t$ 986 # We have $L$ in `transform_cholesky`, so we want to do a linear 987 # solve of $L x = e_t$ where L is lower triangular 988 lapack.{{prefix}}trtrs("L", "N", "N", &self._k_endog, &inc, 989 &self.transform_cholesky[0,0], &self._k_endog, 990 &self.selected_obs[0], &self._k_endog, 991 &info) 992 993 # Calculate loglikelihood contribution of this observation 994 995 # $e_t' H_t^{-1} e_t = (L^{-1} e_t)' L^{-1} e_t = \sum_i e_{i,t}**2$ 996 self.collapse_loglikelihood = 0 997 for i in range(self._k_endog): 998 self.collapse_loglikelihood = self.collapse_loglikelihood + self.selected_obs[i]**2 999 1000 # (p-m) log( 2*pi) + log( |H_t| ) 1001 self.collapse_loglikelihood = ( 1002 self.collapse_loglikelihood + 1003 (self._k_endog - k_states)*{{combined_prefix}}log(2*NPY_PI) + 1004 self.transform_determinant 1005 ) 1006 1007 # -0.5 * ... 1008 self.collapse_loglikelihood = -0.5 * self.collapse_loglikelihood 1009 1010 # Set pointers 1011 self._obs = &self.collapse_obs[0] 1012 self._design = &self.collapse_design[0,0] 1013 self._obs_cov = &self.collapse_obs_cov[0,0] 1014 1015 # TODO can I replace this with k_states? I think I should be able to 1016 return self._k_states 1017 1018# ### Selected covariance matrice 1019cdef int {{prefix}}select_cov(int k, int k_posdef, 1020 {{cython_type}} * tmp, 1021 {{cython_type}} * selection, 1022 {{cython_type}} * cov, 1023 {{cython_type}} * selected_cov): 1024 cdef: 1025 {{cython_type}} alpha = 1.0 1026 {{cython_type}} beta = 0.0 1027 1028 # Only need to do something if there is a covariance matrix 1029 # (i.e k_posdof == 0) 1030 if k_posdef > 0: 1031 1032 # #### Calculate selected state covariance matrix 1033 # $Q_t^* = R_t Q_t R_t'$ 1034 # 1035 # Combine a selection matrix and a covariance matrix to get 1036 # a simplified (but possibly singular) "selected" covariance 1037 # matrix (see e.g. Durbin and Koopman p. 43) 1038 1039 # `tmp0` array used here, dimension $(m \times r)$ 1040 1041 # TODO this does not require two ?gemm calls, since we know that it 1042 # is just selection rows and columns of the Q matrix 1043 1044 # $\\#_0 = 1.0 * R_t Q_t$ 1045 # $(m \times r) = (m \times r) (r \times r)$ 1046 blas.{{prefix}}gemm("N", "N", &k, &k_posdef, &k_posdef, 1047 &alpha, selection, &k, 1048 cov, &k_posdef, 1049 &beta, tmp, &k) 1050 # $Q_t^* = 1.0 * \\#_0 R_t'$ 1051 # $(m \times m) = (m \times r) (m \times r)'$ 1052 blas.{{prefix}}gemm("N", "T", &k, &k, &k_posdef, 1053 &alpha, tmp, &k, 1054 selection, &k, 1055 &beta, selected_cov, &k) 1056 1057{{endfor}}