Created
September 28, 2016 10:59
-
-
Save yogabonito/784bb4e4a36515c04e2a4518e8d348de to your computer and use it in GitHub Desktop.
VAR bugs: 1) lags=0, 2) select_order Raw
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{ | |
"cells": [ | |
{ | |
"cell_type": "code", | |
"execution_count": 1, | |
"metadata": { | |
"collapsed": true | |
}, | |
"outputs": [], | |
"source": [ | |
"import pandas\n", | |
"from statsmodels.tsa.vector_ar.var_model import VAR\n", | |
"import data as dta # statsmodels/datasets/macrodata/data.py\n", | |
"from statsmodels.tsa.base.datetools import dates_from_str\n", | |
"\n", | |
"# prepare data\n", | |
"mdata = dta.load_pandas().data\n", | |
"iidata = dta.load_pandas();\n", | |
"mdata = iidata.data\n", | |
"dates = mdata[['year', 'quarter']].astype(int).astype(str)\n", | |
"quarterly = dates[\"year\"] + \"Q\" + dates[\"quarter\"]\n", | |
"from statsmodels.tsa.base.datetools import dates_from_str\n", | |
"quarterly = dates_from_str(quarterly)\n", | |
"mdata = mdata[[\"realcons\", \"realgdp\", \"realinv\"]]\n", | |
"mdata.index = pandas.DatetimeIndex(quarterly)\n", | |
"data = mdata" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"# Bug 1" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"If `lags > 0`, everything is OK:" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 2, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [], | |
"source": [ | |
"model = VAR(data)\n", | |
"est = model.fit(maxlags=4, trend=\"nc\")" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"If `lags == 0` but there is some trend, everything is OK too:" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 3, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [], | |
"source": [ | |
"model = VAR(data)\n", | |
"est = model.fit(maxlags=0, trend=\"c\")" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"If `lags == 0` and there is no trend, we get an exception. The reason why we get this exception is that the matrix Z (see Lütkepohl, p. 70) has shape (203, 0) in this case, i.e. it doesn't have columns." | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 4, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"ename": "ValueError", | |
"evalue": "math domain error", | |
"output_type": "error", | |
"traceback": [ | |
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", | |
"\u001b[1;31mValueError\u001b[0m Traceback (most recent call last)", | |
"\u001b[1;32m<ipython-input-4-6ffc99743af7>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[0mmodel\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mVAR\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0mest\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mmodel\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfit\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmaxlags\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtrend\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m\"nc\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m", | |
"\u001b[1;32m/home/vegcev/D/TUG/Masterarbeit/GSoC/statsmodels-yogabonito/statsmodels/tsa/vector_ar/var_model.py\u001b[0m in \u001b[0;36mfit\u001b[1;34m(self, maxlags, method, ic, trend, verbose)\u001b[0m\n\u001b[0;32m 436\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mnobs\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mendog\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mlags\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 437\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 438\u001b[1;33m \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_estimate_var\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mlags\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtrend\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mtrend\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 439\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 440\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m_estimate_var\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mlags\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0moffset\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtrend\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'c'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", | |
"\u001b[1;32m/home/vegcev/D/TUG/Masterarbeit/GSoC/statsmodels-yogabonito/statsmodels/tsa/vector_ar/var_model.py\u001b[0m in \u001b[0;36m_estimate_var\u001b[1;34m(self, lags, offset, trend)\u001b[0m\n\u001b[0;32m 459\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 460\u001b[0m \u001b[1;31m# Lutkepohl p75, about 5x faster than stated formula\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 461\u001b[1;33m \u001b[0mparams\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mlstsq\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mz\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0my_sample\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 462\u001b[0m \u001b[0mresid\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0my_sample\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mz\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mparams\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 463\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n", | |
"\u001b[1;32m/home/vegcev/.local/lib/python3.4/site-packages/numpy/linalg/linalg.py\u001b[0m in \u001b[0;36mlstsq\u001b[1;34m(a, b, rcond)\u001b[0m\n\u001b[0;32m 1887\u001b[0m \u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbstar\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0m_to_native_byte_order\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbstar\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 1888\u001b[0m \u001b[0ms\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mzeros\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmin\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mm\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mn\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mreal_t\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1889\u001b[1;33m \u001b[0mnlvl\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mmax\u001b[0m\u001b[1;33m(\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mint\u001b[0m\u001b[1;33m(\u001b[0m \u001b[0mmath\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mlog\u001b[0m\u001b[1;33m(\u001b[0m \u001b[0mfloat\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmin\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mm\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mn\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m/\u001b[0m\u001b[1;36m2.\u001b[0m \u001b[1;33m)\u001b[0m \u001b[1;33m)\u001b[0m \u001b[1;33m+\u001b[0m \u001b[1;36m1\u001b[0m \u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 1890\u001b[0m \u001b[0miwork\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mzeros\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m3\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mmin\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mm\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mn\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mnlvl\u001b[0m\u001b[1;33m+\u001b[0m\u001b[1;36m11\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mmin\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mm\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mn\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfortran_int\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 1891\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0misComplexType\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mt\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", | |
"\u001b[1;31mValueError\u001b[0m: math domain error" | |
] | |
} | |
], | |
"source": [ | |
"model = VAR(data)\n", | |
"est = model.fit(maxlags=0, trend=\"nc\")" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"# Bug 2" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"`select_order()` doesn't take into account any deterministic terms (it always calls `_estimate_var()` without trend argument). So `_estimate_var()` uses the default trend `\"c\"`. If we add a `trend` parameter to `select_order()`, we can pass it to `_estimate_var()`. (In case of `trend == \"nc\"`, this will raise the same exception as above because `select_order()` will call `_estimate_var()` with `lags=0`.)" | |
] | |
} | |
], | |
"metadata": { | |
"kernelspec": { | |
"display_name": "Python 3", | |
"language": "python", | |
"name": "python3" | |
}, | |
"language_info": { | |
"codemirror_mode": { | |
"name": "ipython", | |
"version": 3 | |
}, | |
"file_extension": ".py", | |
"mimetype": "text/x-python", | |
"name": "python", | |
"nbconvert_exporter": "python", | |
"pygments_lexer": "ipython3", | |
"version": "3.4.3+" | |
} | |
}, | |
"nbformat": 4, | |
"nbformat_minor": 0 | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment