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

[LakeModel] Simplify computation of the steady state #355

Merged
merged 11 commits into from
Dec 20, 2023

Conversation

mnshkw
Copy link
Contributor

@mnshkw mnshkw commented Jul 10, 2023

Compute the steady state by the straightforward formula for a 2x2 Markov chain instead of iteration.
See also: #169.

@jstac
Copy link
Contributor

jstac commented Jul 10, 2023

Thanks @mnshkw , appreciated.

@shlff , could you please do a first round review of this and also consider #169 ?

@jstac
Copy link
Contributor

jstac commented Jul 10, 2023

Thanks @mnshkw , appreciated.

@shlff , could you please do a first round review of this and also consider #169 ?

The build error might need the attention of @mmcky

error = np.max(np.abs(new_x - x))
x = new_x
x = np.array([self.A_hat[0, 1], self.A_hat[1, 0]])
x /= x.sum()
Copy link
Member

Choose a reason for hiding this comment

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

Suggest

        x /= x.sum()
        return x

to

        return x / x.sum()

Copy link
Member

Choose a reason for hiding this comment

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

The former is in principle strictly better: x / x.sum() creates an intermediate array while x /= x.sum() is an in-place operation. (Of course, the difference is negligible for this size of the data.)

@shlff
Copy link
Member

shlff commented Jul 11, 2023

Thanks @mnshkw and @jstac .

I quickly reviewed this change and issue #169 , which makes the most of the analytical solution of the stationary distribution of the Markov chain when the stochastic matrix is positive.

For example, let the stochastic matrix be

$$ P = \left( \begin{matrix} 1 - \lambda & \lambda \\ \alpha & 1 - \alpha \end{matrix} \right) $$

if $\alpha \in (0, 1)$ and $\lambda \in (0, 1)$, then $P$ has a unique stationary distribution.

It is a smart change, but the thing is that it cannot handle the case when $\alpha$ or $\lambda$ takes the boundary values, that is, $\alpha, \lambda$ take $0, 1$.

I suggest that we turn this change or the original one into an exercise.

@oyamad
Copy link
Member

oyamad commented Jul 11, 2023

@shlff Good point. The proposed code does not work if lambda = alpha = 0 (in your notation) (it does as long as at least one is strictly positive). But in this case, there is nothing to analyze?

@shlff
Copy link
Member

shlff commented Jul 11, 2023

Thanks for clarifying it @oyamad . Sorry for my late reply.

I reckon there are two benefits of mentioning the extreme case (that is, lambda = alpha = 0).

One is that we can treat the case as a limiting reference point for our analysis of varying values of lambda, alpha in lecture lake_model.

Additionally, mentioning this case enhances the audience's understanding of how we can apply the improved method, when it will fail, and how we can resolve the issue when it fails.

@oyamad
Copy link
Member

oyamad commented Jul 12, 2023

@shlff Thanks, but I am afraid I could not see your point. The arguments in that lecture all rely on the assumptions of irreducibility (lambda >0, alpha > 0) and aperiodicity ((lambda, alpha) neq (1, 1)).

The formula psi^*[0] = alpha / (alpha + lambda) fails if and only if alpha = 0 and lambda = 0. In this case, any distribution is obviously a stationary distribution?

Could you elaborate?

@shlff
Copy link
Member

shlff commented Jul 13, 2023

The following part of the lecture lake_model could be simplified similarly:

def rate_steady_state(self, tol=1e-6):
"""
Finds the steady state of the system :math:`x_{t+1} = \hat A x_{t}`
Returns
--------
xbar : steady state vector of employment and unemployment rates
"""
x = np.full(2, 0.5)
error = tol + 1
while error > tol:
new_x = self.A_hat @ x
error = np.max(np.abs(new_x - x))
x = new_x
return x

Also, consider using inheritance for the class.

@mmcky
Copy link
Contributor

mmcky commented Aug 8, 2023

Execution issue is due to: https://github.com/orgs/QuantEcon/discussions/108

@mmcky
Copy link
Contributor

mmcky commented Aug 9, 2023

I need to fix #359 to fix the build

@mmcky
Copy link
Contributor

mmcky commented Dec 13, 2023

@shlff would you mind to have a final review of this PR? Is anything outstanding?

Preview builds will not run as the source change is from an external fork.

@mmcky mmcky requested a review from shlff December 13, 2023 23:11
@shlff
Copy link
Member

shlff commented Dec 14, 2023

Thanks @mmcky . I will check this PR and get back to you by Friday 15th.

@shlff
Copy link
Member

shlff commented Dec 18, 2023

Thanks @mmcky . After thoroughly thinking of the simplified approach and related lectures, we should adopt this change.

  • Additionally, I suggest that we add a note here to remind readers that we will discuss in the section ergodicity why we can calculate the steady state in this way and the conditions we should check to ensure the existence.
  • Furthermore, we should also apply this change to the corresponding method of a class in one of the exercises of that lecture:

def rate_steady_state(self, tol=1e-6):
"""
Finds the steady state of the system :math:`x_{t+1} = \hat A x_{t}`
Returns
--------
xbar : steady state vector of employment and unemployment rates
"""
x = np.full(2, 0.5)
error = tol + 1
while error > tol:
new_x = self.A_hat @ x
error = np.max(np.abs(new_x - x))
x = new_x
return x

Looking forward to your comments.

@jstac
Copy link
Contributor

jstac commented Dec 18, 2023

Thanks @shlff , please go ahead.

@shlff shlff added the in-work label Dec 19, 2023
@oyamad
Copy link
Member

oyamad commented Dec 19, 2023

I am afraid I did not understand this argument:

Additionally, I suggest that we add a note here to remind readers that we will discuss in the section ergodicity why we can calculate the steady state in this way and the conditions we should check to ensure the existence.

The code comes from the formula "inflow = outflow" at a steady state: x[1] * P[1, 0] = x[0] * P[0, 1], with the sum being 1. It has nothing to do with "ergodicity"; it is just the definition of a stationary distribution.

The reference should be to "25.6. Stationary Distributions" instead.

@jstac
Copy link
Contributor

jstac commented Dec 19, 2023

I agree with @oyamad .

@shlff
Copy link
Member

shlff commented Dec 19, 2023

The code comes from the formula "inflow = outflow" at a steady state: x[1] * P[1, 0] = x[0] * P[0, 1], with the sum being 1. It has nothing to do with "ergodicity"; it is just the definition of a stationary distribution.

The reference should be to "25.6. Stationary Distributions" instead.

Thanks @oyamad . I am sorry that's a typo. Actually, I refer it to 69.4. Dynamics of an Individual Worker, where there is a discussion on stationary distribution.

lectures/lake_model.md Outdated Show resolved Hide resolved
@jstac
Copy link
Contributor

jstac commented Dec 19, 2023

@mmcky , would you mind to fix the build?

@mmcky
Copy link
Contributor

mmcky commented Dec 19, 2023

@jstac this build cannot be fixed as it is originating from a fork (GitHub doesn't want to give out our EC2 credentials -- which is actually nice of them :-)) but it is a pain. I haven't been able to think of a good technical solution to this problem yet. It's only an issue for lectures that make use of EC2 as the backend (i.e. GPU).

I will make a replicate of the PR on a local branch so you can preview.

mmcky added a commit that referenced this pull request Dec 19, 2023
@mmcky mmcky mentioned this pull request Dec 19, 2023
mmcky and others added 2 commits December 19, 2023 16:02
Co-authored-by: Daisuke Oyama <[email protected]>
Co-authored-by: Daisuke Oyama <[email protected]>
lectures/lake_model.md Outdated Show resolved Hide resolved
@oyamad
Copy link
Member

oyamad commented Dec 19, 2023

@shlff Sorry you are right! A_hat is column-stochastic (not row-stochastic as in P in Lecture 25).

So the formula "inflow = outflow" here is: A_hat[0, 1] * x[1] = A_hat[1, 0] * x[0].

@oyamad
Copy link
Member

oyamad commented Dec 19, 2023

2473f00 and d7d8f3a should be reverted, sorry...

@mmcky
Copy link
Contributor

mmcky commented Dec 19, 2023

thanks @oyamad and @shlff

@mmcky mmcky added review and removed in-work labels Dec 19, 2023
@mmcky
Copy link
Contributor

mmcky commented Dec 19, 2023

@mmcky mmcky mentioned this pull request Dec 19, 2023
@mmcky mmcky requested a review from jstac December 19, 2023 22:31
@jstac
Copy link
Contributor

jstac commented Dec 19, 2023

Thanks @oyamad @mmcky @shlff .

  • "steady state of the rate" -> "steady state rate"

  • steady state rate needs to be defined

  • "(which we discuss in section Dynamics of an Individual Worker)" should be "using [a technique](same link) we previously introduced for computing stationary distributions of Markov chains"

@mmcky I'll step out of this now so please merge once these changes are made and you are happy (or request someone else to review).

@mmcky
Copy link
Contributor

mmcky commented Dec 19, 2023

steady state rate needs to be defined

@shlff can I get your help in addressing this comment?

thanks for your review @jstac

@mmcky
Copy link
Contributor

mmcky commented Dec 19, 2023

@shlff
Copy link
Member

shlff commented Dec 19, 2023

steady state rate needs to be defined

@shlff can I get your help in addressing this comment?

thanks for your review @jstac

Thanks @mmcky . I found the steady state rate has already been defined in the texts below the code:

On the other hand, the vector of employment and unemployment rates $x_t$ can be in a steady state $\bar x$ if
there exists an $\bar x$ such that
* $\bar x = \hat A \bar x$
* the components satisfy $\bar e + \bar u = 1$
This equation tells us that a steady state level $\bar x$ is an eigenvector of $\hat A$ associated with a unit eigenvalue.

I suggest we either leave it as it is or shift it up to the end of

we can also write this as
$$
x_{t+1} = \hat A x_t
\quad \text{where} \quad
\hat A := \frac{1}{1 + g} A
$$
You can check that $e_t + u_t = 1$ implies that $e_{t+1}+u_{t+1} = 1$.
This follows from the fact that the columns of $\hat A$ sum to 1.

lectures/lake_model.md Outdated Show resolved Hide resolved
Co-authored-by: mmcky <[email protected]>
lectures/lake_model.md Outdated Show resolved Hide resolved
@mmcky
Copy link
Contributor

mmcky commented Dec 20, 2023

thanks for your comments @shlff, @oyamad, and @jstac. I will merge this now.

Thanks @mnshkw for your PR -- greatly appreciated.

@mmcky mmcky merged commit b617df5 into QuantEcon:main Dec 20, 2023
5 of 6 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants