Merge 24afab2 into 8ddede8

examples/notebooks/statespace_varmax.ipynb

@@ -0,0 +1,140 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# VARMAX models\n",
+    "\n",
+    "This is a notebook stub for VARMAX models. Full development will be done after impulse response functions are available."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import statsmodels.api as sm\n",
+    "import matplotlib.pyplot as plt"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "dta = sm.datasets.webuse('lutkepohl2', 'http://www.stata-press.com/data/r12/')\n",
+    "dta.index = dta.qtr\n",
+    "endog = dta.ix['1960-04-01':'1978-10-01', ['dln_inv', 'dln_inc', 'dln_consump']]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Model specification\n",
+    "\n",
+    "The `VARMAX` class in Statsmodels allows estimation of VAR, VMA, and VARMA models (through the `order` argument), optionally with a constant term (via the `trend` argument). Exogenous regressors may also be included (as usual in Statsmodels, by the `exog` argument), and in this way a time trend may be added. Finally, the class allows measurement error (via the `measurement_error` argument) and allows specifying either a diagonal or unstructured innovation covariance matrix (via the `error_cov_type` argument)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Example 1: VAR\n",
+    "\n",
+    "Below is a simple VARX(2) model in two endogenous variables and an exogenous series, but no constant term. Notice that we needed to allow for more iterations than the default (which is `maxiter=50`) in order for the likelihood estimation to converge. This is not unusual in VAR models which have to estimate a large number of parameters, often on a relatively small number of time series: this model, for example, estimates 27 parameters off of 75 observations of 3 variables."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "# exog = pd.Series(np.arange(len(endog)), index=endog.index, name='trend')\n",
+    "exog = endog['dln_consump']\n",
+    "mod = sm.tsa.VARMAX(endog[['dln_inv', 'dln_inc']], order=(2,0), trend='nc', exog=exog)\n",
+    "res = mod.fit(maxiter=1000)\n",
+    "print res.summary()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Example 2: VMA\n",
+    "\n",
+    "A vector moving average model can also be formulated. Below we show a VMA(2) on the same data, but where the innovations to the process are uncorrelated. In this example we leave out the exogenous regressor but now include the constant term."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "mod = sm.tsa.VARMAX(endog[['dln_inv', 'dln_inc']], order=(0,2), error_cov_type='diagonal')\n",
+    "res = mod.fit(maxiter=1000)\n",
+    "print res.summary()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Caution: VARMA(p,q) specifications\n",
+    "\n",
+    "Although the model allows estimating VARMA(p,q) specifications, these models are not identified without additional restrictions on the representation matrices, which are not built-in. For this reason, it is recommended that the user proceed with error (and indeed a warning is issued when these models are specified). Nonetheless, they may in some circumstances provide useful information."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "mod = sm.tsa.VARMAX(endog[['dln_inv', 'dln_inc']], order=(1,1))\n",
+    "res = mod.fit(maxiter=1000)\n",
+    "print res.summary()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 2",
+   "language": "python",
+   "name": "python2"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 2
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython2",
+   "version": "2.7.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}

statsmodels/src/blas_lapack.pxd

@@ -593,6 +593,23 @@ ctypedef int spotrs_t(
     int *info         # 0 if success, otherwise an error code (integer)
 )
 
+ctypedef int strtrs_t(
+    #  STRTRS solves a triangular system of the form
+    #      A * X = B,  A**T * X = B,  or  A**H * X = B,
+    # where A is a triangular matrix of order N, and B is an N-by-NRHS
+    # matrix.  A check is made to verify that A is nonsingular.
+    char *uplo,          # 'U':  A is upper triangular
+    char *trans,         # N: A * X = B; T: A**T * X = B; C: A**H * X = B
+    char *diag,          # {'U','N'}: unit triangular or not
+    int *n,              # The order of the matrix A.  n >= 0.
+    int *nrhs,           # The number of right hand sides
+    np.float32_t *a,     # Matrix A: nxn
+    int *lda,            # The size of the first dimension of A (in memory)
+    np.float32_t *b,     # Matrix B: nxnrhs
+    int *ldb,            # The size of the first dimension of B (in memory)
+    int *info            # 0 if success, otherwise an error code (integer)
+)
+
 ctypedef int dgetrf_t(
     # DGETRF - compute an LU factorization of a general M-by-N
     # matrix A using partial pivoting with row interchanges
@@ -666,6 +683,23 @@ ctypedef int dpotrs_t(
     int *info         # 0 if success, otherwise an error code (integer)
 )
 
+ctypedef int dtrtrs_t(
+    #  DTRTRS solves a triangular system of the form
+    #      A * X = B,  A**T * X = B,  or  A**H * X = B,
+    # where A is a triangular matrix of order N, and B is an N-by-NRHS
+    # matrix.  A check is made to verify that A is nonsingular.
+    char *uplo,          # 'U':  A is upper triangular
+    char *trans,         # N: A * X = B; T: A**T * X = B; C: A**H * X = B
+    char *diag,          # {'U','N'}: unit triangular or not
+    int *n,              # The order of the matrix A.  n >= 0.
+    int *nrhs,           # The number of right hand sides
+    np.float64_t *a,     # Matrix A: nxn
+    int *lda,            # The size of the first dimension of A (in memory)
+    np.float64_t *b,     # Matrix B: nxnrhs
+    int *ldb,            # The size of the first dimension of B (in memory)
+    int *info            # 0 if success, otherwise an error code (integer)
+)
+
 ctypedef int cgetrf_t(
     # CGETRF - compute an LU factorization of a general M-by-N
     # matrix A using partial pivoting with row interchanges
@@ -739,6 +773,23 @@ ctypedef int cpotrs_t(
     int *info           # 0 if success, otherwise an error code (integer)
 )
 
+ctypedef int ctrtrs_t(
+    #  CTRTRS solves a triangular system of the form
+    #      A * X = B,  A**T * X = B,  or  A**H * X = B,
+    # where A is a triangular matrix of order N, and B is an N-by-NRHS
+    # matrix.  A check is made to verify that A is nonsingular.
+    char *uplo,          # 'U':  A is upper triangular
+    char *trans,          # N: A * X = B; T: A**T * X = B; C: A**H * X = B
+    char *diag,           # {'U','N'}: unit triangular or not
+    int *n,              # The order of the matrix A.  n >= 0.
+    int *nrhs,           # The number of right hand sides
+    np.complex64_t *a,   # Matrix A: nxn
+    int *lda,            # The size of the first dimension of A (in memory)
+    np.complex64_t *b,   # Matrix B: nxnrhs
+    int *ldb,            # The size of the first dimension of B (in memory)
+    int *info            # 0 if success, otherwise an error code (integer)
+)
+
 ctypedef int zgetrf_t(
     # ZGETRF - compute an LU factorization of a general M-by-N
     # matrix A using partial pivoting with row interchanges
@@ -810,4 +861,21 @@ ctypedef int zpotrs_t(
     np.complex128_t *b,  # Matrix B: nxnrhs
     int *ldb,            # The size of the first dimension of B (in memory)
     int *info            # 0 if success, otherwise an error code (integer)
+)
+
+ctypedef int ztrtrs_t(
+    #  ZTRTRS solves a triangular system of the form
+    #      A * X = B,  A**T * X = B,  or  A**H * X = B,
+    # where A is a triangular matrix of order N, and B is an N-by-NRHS
+    # matrix.  A check is made to verify that A is nonsingular.
+    char *uplo,          # 'U':  A is upper triangular
+    char *trans,          # N: A * X = B; T: A**T * X = B; C: A**H * X = B
+    char *diag,           # {'U','N'}: unit triangular or not
+    int *n,              # The order of the matrix A.  n >= 0.
+    int *nrhs,           # The number of right hand sides
+    np.complex128_t *a,  # Matrix A: nxn
+    int *lda,            # The size of the first dimension of A (in memory)
+    np.complex128_t *b,  # Matrix B: nxnrhs
+    int *ldb,            # The size of the first dimension of B (in memory)
+    int *info            # 0 if success, otherwise an error code (integer)
 )

statsmodels/tsa/api.py

@@ -17,4 +17,6 @@
 from .x13 import x13_arima_select_order
 from .x13 import x13_arima_analysis
 from .statespace import api as statespace
-from .statespace.sarimax import SARIMAX
+from .statespace.sarimax import SARIMAX
+from .statespace.varmax import VARMAX
+from .statespace.dynamic_factor import DynamicFactor