Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Extend LinearStateSpace class #569

Merged
merged 4 commits into from
Apr 16, 2021
Merged

Extend LinearStateSpace class #569

merged 4 commits into from
Apr 16, 2021

Conversation

shizejin
Copy link
Member

Modify the LinearStateSpace class with two new features:

  1. to compute the steady state covariance matrix by solving a discrete Lyapunov equation rather than "iterating to convergence"

    • the stationary_distributions now does not take arguments max_iter nor tol
    • when I was trying to apply this new method to some examples, I found that the doubling algorithm does not always work (see more details in Doubling algorithm not working for solving discrete Lyapunov equation #568 ]). Therefore, I am forcing stationary_distributions to use the Bartels-Stewart algorithm.
  2. adding a method that will compute the stationary covariance matrix
    E (y_t - \mu_y) (x_t -\mu_x)' = G \Sigma_x(\infty).

    • this is returned as the last output of stationary_distributions method
    • it is also stored as an field Sigma_yx

@pep8speaks
Copy link

pep8speaks commented Jan 16, 2021

Hello @shizejin! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

Line 451:5: E125 continuation line with same indent as next logical line
Line 451:13: E101 indentation contains mixed spaces and tabs
Line 451:13: W191 indentation contains tabs

Line 55:80: E501 line too long (95 > 79 characters)

Comment last updated at 2021-02-21 16:09:12 UTC

@coveralls
Copy link

coveralls commented Jan 16, 2021

Coverage Status

Coverage decreased (-0.03%) to 94.133% when pulling f234eea on shizejin:ext_lss into 02a63b3 on QuantEcon:master.

Copy link
Member

@QBatista QBatista left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@shizejin Thanks for the PR! I'll check if it works for the Townsend lecture.


return mu_x_star, mu_y_star, Sigma_x_star, Sigma_y_star
return mu_x, mu_y, Sigma_x, Sigma_y, Sigma_yx
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This will probably break a significant amount of code using this function. An easy way around this is to use a keyword argument to optionally return Sigma_yx. Otherwise, I'm guessing we'll need a PR to update multiple lectures.

r"""
Reorder the states by shifting the constant terms to the
top of the state vector. Then partition the linear state
space model following Appendix C in RMT Ch2 such that the
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this Appendix C in the latest version of RMT? In an older version, it seems that it was "Appendix B". If so, we might need a more precise reference.

quantecon/lss.py Outdated
A_diag = np.diag(A)
num_const = 0
for idx in range(n):
if (A_diag[idx] == 1) and (C[idx, :] == 0).all():
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This doesn't seem to be a sufficient condition for finding constants. For instance, consider:

A = np.array([[1., -0.1], [0.5, 0.9]])
C = np.array([[0.], [1.]])

This has no constants and A is apparently stable since:

np.absolute(np.linalg.eig(A)[0])

yields

array([0.97467943, 0.97467943])

Do we really need re-ordering, or can we ask the user to provide a system that takes the special form 2.C.1?

Detecting constants in general is tricky. Consider:

A = np.array([[0.5, 0.5], [0.5, 0.5]])
C = np.array([[0.], [0.]])

If the initial conditions are the same, then this produces constant processes.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree that this part is actually tricky. But I guess we need to do this inside the function so that we don't complicate the input format, with the same logic you mentioned above.

Essentially we want to detect those states associated with unity eigenvalues (we do not allow for eigenvalues strictly greater than 1). In addition to the current two conditions for finding the targeted states, I guess we can check np.linalg.norm(row_vec)==1 or we can compute the eigenvalues directly.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The second example I gave has eigenvalues 1 and 0, and none of the rows satisfy np.linalg.norm(row_vec)==1. Besides, a row may satisfy np.linalg.norm(row_vec)==1, but not be a vector containing one 1 and 0s otherwise.

I think the following procedure might be good enough, even if it doesn't cover all possible cases:

  1. Check whether the ith row is a vector with a 1 in the ith entry, and 0s otherwise. Store these rows as "constant" rows.
  2. Re-order the matrices based on 1.
  3. Verify that A_22 is stable by computing the modulus of eigenvalues. If the check fails, raise an error.

What do you think?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Isn't step 1 equivalent to (A_diag[idx] == 1) and np.linalg.norm(row_vec)==1? Sorry that I might not clear in my previous reply, I meant to add np.linalg.norm(row_vec)==1 "additionally". And (C[idx, :] == 0).all() is also a necessary condition to check for a constant state.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, this looks good to me. I thought you meant only checking np.linalg.norm(row_vec)==1.

@QBatista
Copy link
Member

This method works for the Townsend lecture -- thanks!

@mmcky
Copy link
Contributor

mmcky commented Jan 17, 2021

Thanks @shizejin. @QBatista one you're happy with this PR would you be willing to approve?

@shizejin we will need to think of backward compatibility. Is this just two new methods (i.e. will it impact existing users of the function i.e. through an altered interface)? It seems to be modifying existing method interfaces etc.

@mmcky mmcky marked this pull request as draft January 17, 2021 23:46
@QBatista
Copy link
Member

@mmcky Yes, sure.

@shizejin
Copy link
Member Author

Thanks @shizejin. @QBatista one you're happy with this PR would you be willing to approve?

@shizejin we will need to think of backward compatibility. Is this just two new methods (i.e. will it impact existing users of the function i.e. through an altered interface)? It seems to be modifying existing method interfaces etc.

Yes @mmcky this PR is altering the existing method stationary_distributions. Especially, I added one more stationary moment to the output of this method, and two keyword argument inputs max_iter and tol are removed. The first change is more drastic because it means all the existing code that calls this method by mu_x, mu_y, Sigma_x, Sigma_y = lss.stationary_distributions() will fail due to the extra output Sigma_yx to unpack. But given that we want to have this moment returned by stationary_distributions, this is the most appropriate way to modify the existing method. Would you think it to be feasible to adjust all the lectures which calls this lss.stationary_distributions to adopt this modification?

@jstac
Copy link
Contributor

jstac commented Jan 19, 2021

@shizejin It's not a big deal to change the lecture code to handle the extra output, although I do wonder if there are people using this routine in the wild who might be inconvenienced...?

There are various alternatives, such as a flag like compute_covariance=True when the extra output is desired, but I don't particularly like them.

Overall, this seems to me like a valuable enhancement and it would be best to just change the lecture code to accommodate it. I'm happy to do that @mmcky and @oyamad .

@mmcky
Copy link
Contributor

mmcky commented Jan 19, 2021

thanks @jstac -- I think breaking changes are reasonable until we issue 1.0 but @shizejin can you please document this change clearly with a deprecation notice (i..e from and to) so that I can make this clear in the next release notes.

@shizejin
Copy link
Member Author

thanks @jstac -- I think breaking changes are reasonable until we issue 1.0 but @shizejin can you please document this change clearly with a deprecation notice (i..e from and to) so that I can make this clear in the next release notes.

Thanks @mmcky and @jstac for very useful comments and suggestions. I will write down a clear note once I modify this PR adopting some suggestions from @QBatista . Thank you all!

Add conditions for finding the constant term in the transition matrix.
@shizejin
Copy link
Member Author

I have modified the code according to the great suggestion by @QBatista.

@mmcky would you please let me know if the following description looks clear enough to you?

The two keyword arguments, tol and max_iter, of qe.LinearStateSpace.stationary_distributions have been deprecated. Now it does NOT take any keyword argument as we do not use the iteration method to find the stationary values. Instead, we solve a discrete Lyapunov equation.

We also make qe.LinearStateSpace.stationary_distributions to return one more useful output Sigma_yx, which is the stationary covariance matrix between y_x and x_t.

Therefore, the recommended use of qe.LinearStateSpace.stationary_distributions is updated from

mu_x, mu_y, Sigma_x, Sigma_y = lss.stationary_distributions(max_iter, tol)

to

mu_x, mu_y, Sigma_x, Sigma_y, Sigma_yx = lss.stationary_distributions()

Please let me know if anything is unclear to you and if there is any information I could provide :)

@QBatista
Copy link
Member

The code looks good to me. Thanks @shizejin !

@mmcky
Copy link
Contributor

mmcky commented Feb 23, 2021

thanks @shizejin and @QBatista this looks great.
We will need to put together coordinated updates for code on all quantecon lectures prior to release via PyPI.

@mmcky mmcky marked this pull request as ready for review February 23, 2021 03:43
@mmcky
Copy link
Contributor

mmcky commented Mar 26, 2021

@shizejin @QBatista will you coordinate PR's for the lectures? Or do I need to request some help from @shlff?

@shizejin
Copy link
Member Author

Hi @mmcky, this should be easy. I will try to wrap up a PR by this weekend. If I encounter some difficulties, I will reach @QBatista and @shlff for help. What do you think?

@mmcky
Copy link
Contributor

mmcky commented Mar 26, 2021

Thanks @shizejin that would be great.

@shizejin
Copy link
Member Author

Hi @mmcky could you please let me know if the two PRs look good to you?

@mmcky
Copy link
Contributor

mmcky commented Apr 8, 2021

Thanks @shizejin -- looking good. I am holding off on merging this as we are in the middle of a migration for:

python.quantecon.org

to the new jupyter book system. Sorry for the delay. I will merge this PR next week and then issue a new release along with the additions to the lecture repos.

@mmcky
Copy link
Contributor

mmcky commented Apr 16, 2021

I am going to release this now -- and then migrate the advanced python lectures PR as we migrate those to jupyter-book

@mmcky
Copy link
Contributor

mmcky commented Apr 16, 2021

thanks @shizejin

@mmcky mmcky merged commit d975adb into QuantEcon:master Apr 16, 2021
mmcky added a commit to QuantEcon/lecture-python.myst that referenced this pull request Apr 28, 2021
* update QuantEcon.py QuantEcon/QuantEcon.py#569

* catch execution errors for CI

* TEST: test against new release of QuantEcon

* add .git for install from quantecon.py@new-release

* use correct branch name

* update environment
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

6 participants