ENH: add rsquared for statespace model #4734 #6620

BenjaminLiuPenrose · 2020-04-04T22:13:20Z

closes ENH: Add rsquared and variants to state space models #4734
tests added / passed.
code/documentation is well formatted.
properly formatted commit message. See
NumPy's guide.

Notes:

It is essential that you add a test when making code changes. Tests are not
needed for doc changes.
When adding a new function, test values should usually be verified in another package (e.g., R/SAS/Stata).
When fixing a bug, you must add a test that would produce the bug in master and
then show that it is fixed with the new code.
New code additions must be well formatted. Changes should pass flake8. If on Linux or OSX, you can
verify you changes are well formatted by running
```
git diff upstream/master -u -- "*.py" | flake8 --diff --isolated
```
assuming flake8 is installed. This command is also available on Windows
using the Windows System for Linux once flake8 is installed in the
local Linux environment. While passing this test is not required, it is good practice and it help
improve code quality in statsmodels.
Docstring additions must render correctly, including escapes and LaTeX.

pep8speaks · 2020-04-04T22:13:24Z

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

There are currently no PEP 8 issues detected in this Pull Request. Cheers! 🍻

Comment last updated at 2020-04-17 23:11:02 UTC

coveralls · 2020-04-04T23:16:34Z

Coverage increased (+0.02%) to 87.707% when pulling 4bead54 on BenjaminLiuPenrose:my_change into 35b04e5 on statsmodels:master.

codecov · 2020-04-04T23:20:38Z

Codecov Report

Merging #6620 into master will decrease coverage by <.01%.
The diff coverage is 25%.

@@            Coverage Diff             @@
##           master    #6620      +/-   ##
==========================================
- Coverage   85.31%    85.3%   -0.01%     
==========================================
  Files         646      646              
  Lines      103922   103926       +4     
  Branches    11311    11311              
==========================================
+ Hits        88657    88658       +1     
- Misses      12806    12809       +3     
  Partials     2459     2459

Impacted Files	Coverage Δ
statsmodels/tsa/statespace/mlemodel.py	`87.14% <25%> (-0.2%)`	⬇️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 35b04e5...7c37307. Read the comment docs.

BenjaminLiuPenrose · 2020-04-05T02:20:03Z

where should I write the test case to pass the codecov/patch

ChadFulton · 2020-04-06T00:03:25Z

Thanks very much for submitting this PR, @BenjaminLiuPenrose! This will be a really nice feature to have.

I have a couple of suggestions:

For computing the SSE, you're currently using the large-sample approximation of (5.5.3), because self.resid corresponds to v_t and not to \tilde v_t (which is self.standardized_forecasts_error). Instead, why don't we use the finite prediction error variance defined in (5.5.2), so that you would compute the SSE as in (5.5.13):
```
d = np.maximum(self.loglikelihood_burn, self.nobs_diffuse)
srss = np.sum(self.standardized_forecasts_error[0, d:]**2)
f_T = self.forecasts_error_cov[0, 0, -1]
sse = f_T * srss
```
Since there are multiple ways to define the R^2 depending on the baseline model used as a comparison, how about defining one new method, get_rsquared that accepts one argument baseline. baseline could accept either "mean", which is what you currently have as the rsquared attribute, but could also accept "rwdrift", which is what you currently have as rsquared_difference. I think the default value should be "rwdrift". Eventually (or in this PR if you like), we could add a seasonal argument to allow specifying the R^2 that Harvey has in (5.5.17).
Then the two properties that you added would just call that method with the argument baseline="mean", and baseline="rwdrift". I also suggest that you rename the attributes in MLEResults to be rsquared_mean, rsquared_rwdrift.
I would suggest having the rsquared property raise a NotImplementedError, with a message letting users know that in state space models there is not a single comparison model that will always be useful, and suggesting that they use rsquared_rwdrift. You can also mention that rsquared_mean exists, but is not generally recommended. I think this will be helpful, since the basic rsquared_mean is not usually so helpful in time series contexts.
We need to make sure that these don't return erroneous numbers for multivariate models. For this PR, if you want to just support univariate models, that's okay. In that case, you should raise a NotImplementedError if self.model.k_endog > 1. If you want to try to support multivariate models, that's great too!

Thanks again!

BenjaminLiuPenrose · 2020-04-06T22:10:39Z

Thanks very much for submitting this PR, @BenjaminLiuPenrose! This will be a really nice feature to have.

I have a couple of suggestions:
For computing the SSE, you're currently using the large-sample approximation of (5.5.3), because self.resid corresponds to v_t and not to \tilde v_t (which is self.standardized_forecasts_error). Instead, why don't we use the finite prediction error variance defined in (5.5.2), so that you would compute the SSE as in (5.5.13):
d = np.maximum(self.loglikelihood_burn, self.nobs_diffuse)
srss = np.sum(self.standardized_forecasts_error[0, d:]**2)
f_T = self.forecasts_error_cov[0, 0, -1] * srss
sse = f_T * rss
Since there are multiple ways to define the R^2 depending on the baseline model used as a comparison, how about defining one new method, get_rsquared that accepts one argument baseline. baseline could accept either "mean", which is what you currently have as the rsquared attribute, but could also accept "rwdrift", which is what you currently have as rsquared_difference. I think the default value should be "rwdrift". Eventually (or in this PR if you like), we could add a seasonal argument to allow specifying the R^2 that Harvey has in (5.5.17).

Then the two properties that you added would just call that method with the argument baseline="mean", and baseline="rwdrift". I also suggest that you rename the attributes in MLEResults to be rsquared_mean, rsquared_rwdrift.

I would suggest having the rsquared property raise a NotImplementedError, with a message letting users know that in state space models there is not a single comparison model that will always be useful, and suggesting that they use rsquared_rwdrift. You can also mention that rsquared_mean exists, but is not generally recommended. I think this will be helpful, since the basic rsquared_mean is not usually so helpful in time series contexts.

We need to make sure that these don't return erroneous numbers for multivariate models. For this PR, if you want to just support univariate models, that's okay. In that case, you should raise a NotImplementedError if self.model.k_endog > 1. If you want to try to support multivariate models, that's great too!
Thanks again!

Sure. Will work on that

ChadFulton · 2020-04-07T13:13:17Z

Note: there was typo in my code computing the sse, but I have corrected it above.

BenjaminLiuPenrose · 2020-04-11T01:24:17Z

@ChadFulton how about now?

ChadFulton · 2020-04-11T18:29:26Z

That's great, thanks! I will have a change to make a couple of minor comments on the code itself hopefully later today. The only major thing remaining is coming up with some unit tests.

BenjaminLiuPenrose · 2020-04-12T04:02:44Z

That's great, thanks! I will have a change to make a couple of minor comments on the code itself hopefully later today. The only major thing remaining is coming up with some unit tests.

sure;)

BenjaminLiuPenrose · 2020-04-13T21:09:26Z

@ChadFulton any updates?

statsmodels/tsa/statespace/mlemodel.py

ChadFulton · 2020-04-14T03:20:16Z

I have a couple of comments, but the main thing now is to get some unit tests. Are you familiar with unit testing?

Do you know of any resources we could test this against, or have you thought about other ways that we could write some tests?

BenjaminLiuPenrose · 2020-04-14T18:24:55Z

I have a couple of comments, but the main thing now is to get some unit tests. Are you familiar with unit testing?

Do you know of any resources we could test this against, or have you thought about other ways that we could write some tests?

@ChadFulton guess I will add test_summary_rsquared in \statsmodels\statsmodels\tsa\statespace\tests\test_mlemodel.py,

but what I can do for now is to check r2 is displayed but I don't know the 'ground true' value for each of r2 wrt to the dummy_model

also, can you also point me out how to extend r2 metric for multi-endog variables problem?

BenjaminLiuPenrose · 2020-04-15T19:24:23Z

@ChadFulton test case added, it is a trivial case

BenjaminLiuPenrose · 2020-04-16T15:57:38Z

@ChadFulton not sure how to print R2 for multivariate case

ChadFulton · 2020-04-17T03:28:16Z

This is looking good, thanks!

not sure how to print R2 for multivariate case

I agree that this is a tough call. Currently for other output, like the test statistics, I just print the list, but that is not very pretty. One option would be to create a table in the summary output that displays the R2 for each endog variable in the multivariate case.

statsmodels/tsa/statespace/mlemodel.py

BenjaminLiuPenrose · 2020-04-19T02:12:39Z

I'm now printing it as np array

BenjaminLiuPenrose · 2020-04-20T22:12:07Z

any ideas why the ci/appveyor/pr test failed?

ChadFulton · 2020-04-21T00:16:15Z

any ideas why the ci/appveyor/pr test failed?

I'm not sure, but it seems like it's unrelated.

On a side note - this is looking good, and I'm sorry for the delay - it will probably take me a day or two to get back to this, based on other time commitments.

BenjaminLiuPenrose · 2020-04-21T21:33:30Z

any ideas why the ci/appveyor/pr test failed?

I'm not sure, but it seems like it's unrelated.

On a side note - this is looking good, and I'm sorry for the delay - it will probably take me a day or two to get back to this, based on other time commitments.

sure, np

BenjaminLiuPenrose · 2020-04-27T01:56:29Z

@ChadFulton updates?

ChadFulton

I've added some more comments. The only big thing that we need to do before this can be merged is to add unit tests for the rwdrift and seasonal cases.

ChadFulton · 2020-04-28T03:32:11Z

statsmodels/tsa/statespace/mlemodel.py

@@ -2894,6 +2894,76 @@ def zvalues(self):
        """
        return self.params / self.bse

+    def get_rsquared(self, baseline="rwdrift", **kwargs):


We can add kwargs later if necessary, but we should avoid it unless we actually need to capture unknown keyword arguments (it can lead to problems with, e.g., misspelled keyword arguments).

statsmodels/tsa/statespace/mlemodel.py

ChadFulton · 2020-04-28T03:42:06Z

statsmodels/tsa/statespace/tests/test_mlemodel.py

+    endog = np.array([1, 2, 4, 8, 16])
+    exog = np.array([1, 2, 4, 8, 16])
+
+    mod = sarimax.SARIMAX(endog, exog, order=(0, 0, 0), trend='c')


I don't think that trend='c' is doing anything here, so we may as well remove it.

statsmodels/tsa/statespace/tests/test_mlemodel.py

ChadFulton · 2020-04-28T03:44:02Z

statsmodels/tsa/statespace/tests/test_mlemodel.py

+    exog = np.array([1, 2, 4, 8, 16])
+
+    mod = sarimax.SARIMAX(endog, exog, order=(0, 0, 0), trend='c')
+    res = mod.fit(disp=-1)


You don't need to fit the model - instead, you could put this after you fit the benchmark model and use:

res = mod.smooth(benchmark_res.params)

ChadFulton · 2020-04-28T04:05:38Z

statsmodels/tsa/statespace/tests/test_mlemodel.py

+def test_summary_rsquared():
+    from statsmodels.regression.linear_model import OLS
+    endog = np.array([1, 2, 4, 8, 16])
+    exog = np.array([1, 2, 4, 8, 16])


This shouldn't be a perfect fitting model, so you should make exog not identical to endog. Also, so that the OLS R^2 matches, exog needs to include a constant column.

ChadFulton · 2020-04-28T04:06:09Z

statsmodels/tsa/statespace/tests/test_mlemodel.py

+    endog = np.array([1, 2, 4, 8, 16])
+    exog = np.array([1, 2, 4, 8, 16])
+
+    mod = sarimax.SARIMAX(endog, exog, order=(0, 0, 0), trend='c')


Add concentrate_scale=True so you don't have to estimate the variance parameter.

BenjaminLiuPenrose · 2020-04-28T17:25:41Z

cases

Yea, I agreed. What do you think will be a good test case for rwdrift and seasonal?

bashtage · 2020-05-14T07:08:09Z

statsmodels/tsa/statespace/mlemodel.py

+            rsquared = 1. - sse / ssm
+        elif baseline == "seasonal":
+            from statsmodels.regression.linear_model import OLS
+            from statsmodels.tools.tools import add_constant


You could just use AutoReg which directly supports this model specification at a high level.

statsmodels.tsa.ar_model.AutoReg(endog, 0, trend='c', seasonal=True, periods=seasonal)

bashtage · 2020-05-14T07:09:33Z

statsmodels/tsa/statespace/mlemodel.py

+        However, we have implmented `rsquared_rwdrift` and `rsquared_mean`
+        It is recommended to use `rsquared_rwdrift`
+        """
+        return self.get_rsquared('')


Doesn't this produce an error?

If the error is intentional, you should directly raise NotImplementedError here. You can pull the common message outside the class to share it.

bashtage · 2020-05-14T07:11:10Z

statsmodels/tsa/statespace/mlemodel.py

+    @cache_readonly
+    def rsquared_mean(self):
+        """
+        (float or array) conventional R-squared, 1 - sse/ssm


The docstring all have the wrong format. Could you please update to NumPy style?

bashtage · 2020-05-14T07:12:50Z

Overall pretty close and a useful contribution. It would be good to get this across the line.

add rsquared for statespace model

8d9e7eb

BenjaminLiuPenrose added 2 commits April 4, 2020 18:15

update tsa.statespace.mlemodel

d9a791e

undo

50fab42

update MLEmodel

87b3e98

BenjaminLiuPenrose added 3 commits April 4, 2020 19:37

update :

fcee4a2

fix patch test failure

c222495

update two methods

7c37307

BenjaminLiuPenrose added 5 commits April 10, 2020 13:39

upgrade r2 calc based on conversation

06e7686

print r2 in mleresult.summary

f1f6547

print r2

d93abaf

typo

31885e0

format

db904ec

BenjaminLiuPenrose added 3 commits April 12, 2020 15:18

add seasonal r2

7b4b6c0

formatting

6a75760

formatting

6fc9640

ChadFulton reviewed Apr 14, 2020

View reviewed changes

statsmodels/tsa/statespace/mlemodel.py Outdated Show resolved Hide resolved

statsmodels/tsa/statespace/mlemodel.py Outdated Show resolved Hide resolved

statsmodels/tsa/statespace/mlemodel.py Outdated Show resolved Hide resolved

reorg code

9b2f8e4

BenjaminLiuPenrose added 5 commits April 15, 2020 15:58

add assert error

0182b65

add assertraise

68c175c

raises

0009afa

assert errors

7147561

adjust line length

7156e4a

ChadFulton reviewed Apr 17, 2020

View reviewed changes

statsmodels/tsa/statespace/mlemodel.py Outdated Show resolved Hide resolved

ChadFulton reviewed Apr 17, 2020

View reviewed changes

statsmodels/tsa/statespace/mlemodel.py Outdated Show resolved Hide resolved

BenjaminLiuPenrose added 4 commits April 17, 2020 12:35

print array

40ad96f

print array

53ddbfe

modify test

f914a87

remove blank lines

4bead54

ChadFulton reviewed Apr 28, 2020

View reviewed changes

ChadFulton added comp-tsa type-enh comp-tsa-statespace and removed comp-tsa labels Apr 28, 2020

bashtage reviewed May 14, 2020

View reviewed changes

ChadFulton mentioned this pull request May 18, 2021

BUG/REF: FRED MD/QD links don't work ChadFulton/tsa-notebooks#26

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

ENH: add rsquared for statespace model #4734 #6620

ENH: add rsquared for statespace model #4734 #6620

BenjaminLiuPenrose commented Apr 4, 2020

pep8speaks commented Apr 4, 2020 •

edited

coveralls commented Apr 4, 2020 •

edited

codecov bot commented Apr 4, 2020 •

edited

BenjaminLiuPenrose commented Apr 5, 2020

ChadFulton commented Apr 6, 2020 •

edited

BenjaminLiuPenrose commented Apr 6, 2020

ChadFulton commented Apr 7, 2020

BenjaminLiuPenrose commented Apr 11, 2020

ChadFulton commented Apr 11, 2020

BenjaminLiuPenrose commented Apr 12, 2020

BenjaminLiuPenrose commented Apr 13, 2020

ChadFulton commented Apr 14, 2020

BenjaminLiuPenrose commented Apr 14, 2020 •

edited

BenjaminLiuPenrose commented Apr 15, 2020

BenjaminLiuPenrose commented Apr 16, 2020

ChadFulton commented Apr 17, 2020

BenjaminLiuPenrose commented Apr 19, 2020

BenjaminLiuPenrose commented Apr 20, 2020

ChadFulton commented Apr 21, 2020

BenjaminLiuPenrose commented Apr 21, 2020

BenjaminLiuPenrose commented Apr 27, 2020

ChadFulton left a comment

ChadFulton Apr 28, 2020

ChadFulton Apr 28, 2020

ChadFulton Apr 28, 2020

ChadFulton Apr 28, 2020

ChadFulton Apr 28, 2020

BenjaminLiuPenrose commented Apr 28, 2020

bashtage May 14, 2020

bashtage May 14, 2020

bashtage May 14, 2020

bashtage May 14, 2020

bashtage commented May 14, 2020

ENH: add rsquared for statespace model #4734 #6620

Are you sure you want to change the base?

ENH: add rsquared for statespace model #4734 #6620

Conversation

BenjaminLiuPenrose commented Apr 4, 2020

pep8speaks commented Apr 4, 2020 • edited

Comment last updated at 2020-04-17 23:11:02 UTC

coveralls commented Apr 4, 2020 • edited

codecov bot commented Apr 4, 2020 • edited

Codecov Report

BenjaminLiuPenrose commented Apr 5, 2020

ChadFulton commented Apr 6, 2020 • edited

BenjaminLiuPenrose commented Apr 6, 2020

ChadFulton commented Apr 7, 2020

BenjaminLiuPenrose commented Apr 11, 2020

ChadFulton commented Apr 11, 2020

BenjaminLiuPenrose commented Apr 12, 2020

BenjaminLiuPenrose commented Apr 13, 2020

ChadFulton commented Apr 14, 2020

BenjaminLiuPenrose commented Apr 14, 2020 • edited

BenjaminLiuPenrose commented Apr 15, 2020

BenjaminLiuPenrose commented Apr 16, 2020

ChadFulton commented Apr 17, 2020

BenjaminLiuPenrose commented Apr 19, 2020

BenjaminLiuPenrose commented Apr 20, 2020

ChadFulton commented Apr 21, 2020

BenjaminLiuPenrose commented Apr 21, 2020

BenjaminLiuPenrose commented Apr 27, 2020

ChadFulton left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

BenjaminLiuPenrose commented Apr 28, 2020

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

bashtage commented May 14, 2020

pep8speaks commented Apr 4, 2020 •

edited

coveralls commented Apr 4, 2020 •

edited

codecov bot commented Apr 4, 2020 •

edited

ChadFulton commented Apr 6, 2020 •

edited

BenjaminLiuPenrose commented Apr 14, 2020 •

edited