Hi @witherscp thanks for opening up an issue!

Checking matrix rank would not be sufficient as there is a linear algebra result "almost every matrix is diagonalizable" (i.e. full rank) because of numerical precision. A better value to check would be the condition number (using numpy/scipy).

Re cross-validation step: I am okay w/ adding this in, but a bunch of open questions that are not that simple. However, what would the score be measured on? Is it how well it fits the training data?, or how well it fits the "next epoch", or "the next n time steps"?

If I check condition number using np.linalg.cond(), how can I know whether it would be considered too large? In my sample data, when there is no common average reference the condition values are around ~1000 for each epoch, whereas they are 10^15 when the common average reference is applied which is ill-conditioned.

In terms of the cross-validation, I was planning on using 5-fold cross validation with: reg = RidgeCV(alphas=np.logspace(-15,5,21), cv=5).fit(z, y_sample) within the function _estimate_var. In this case, I believe the epoch is divided into 5 groups, with each one being left out and the fit determined separately. The score would be based on how well it fits the training data within the epoch. Does that make sense to use? Empirically, in the data I have tried it on, the alpha value tends to be 1e-8 or 1e-7 and the R^2 score within an epoch for the optimal alpha is around 0.4-0.7.

If I check condition number using np.linalg.cond(), how can I know whether it would be considered too large? In my sample data, when there is no common average reference the condition values are around ~1000 for each epoch, whereas they are 10^15 when the common average reference is applied which is ill-conditioned.

We could just come up with a sensible default (e.g. 1e6 or 1e7)?

In terms of the cross-validation, I was planning on using 5-fold cross validation with: reg = RidgeCV(alphas=np.logspace(-15,5,21), cv=5).fit(z, y_sample) within the function _estimate_var. In this case, I believe the epoch is divided into 5 groups, with each one being left out and the fit determined separately. The score would be based on how well it fits the training data within the epoch. Does that make sense to use? Empirically, in the data I have tried it on, the alpha value tends to be 1e-8 or 1e-7 and the R^2 score within an epoch for the optimal alpha is around 0.4-0.7.

Sure I think this makes sense. But I suspect might need to split the epochs such that they predict some validation epoch (might be shorter than the window lengths). But maybe let's try it out in your PR and if you can show a good example on simulated (or real data) that would be good to inform the settings.

To make ill conditioned data in simulation you can just add redundant sensors with small noise or something.

Feel free to open the PR early to get some initial feedback.

Submitted the PR. The dynamic model seems to work well. The residuals are quite small, even less than the OLS. I'm not sure if the average epochs model works as well because the optimal alpha value is the absolute minimum of the np.logspace I provide and if I give a bigger range then the alpha value continues to drop but there is an ill-conditioned solving error. So far, I have only tested with the ECOG dataset using a common average reference.