summary impact measures for spatial autoregressive interaction model using Python

impacts_annotated.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# This notebook provides Python code to compute the summary measure effects for spatial autoregressive interaction model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import pysal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Prepare the scenario from LeSage & Thomas-Agnan (2015)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#this prepares spatial weight for a lattice of 8 units in a single row\n",
    "w = pysal.lat2W(1, 8)\n",
    "w.transform = 'r'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.   1.   0.   0.   0.   0.   0.   0. ]\n",
      " [ 0.5  0.   0.5  0.   0.   0.   0.   0. ]\n",
      " [ 0.   0.5  0.   0.5  0.   0.   0.   0. ]\n",
      " [ 0.   0.   0.5  0.   0.5  0.   0.   0. ]\n",
      " [ 0.   0.   0.   0.5  0.   0.5  0.   0. ]\n",
      " [ 0.   0.   0.   0.   0.5  0.   0.5  0. ]\n",
      " [ 0.   0.   0.   0.   0.   0.5  0.   0.5]\n",
      " [ 0.   0.   0.   0.   0.   0.   1.   0. ]]\n"
     ]
    }
   ],
   "source": [
    "#this prints the underlying row-standardized matrix for the spatial weight object\n",
    "w = w.full()[0]\n",
    "print w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### It can be seen above that the weight matrix prepared here matches that from the publication (p. 195)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Set number of locations (n) and number of flow observations (65)\n",
    "n=8\n",
    "N=64\n",
    "\n",
    "#create n x n identity matrix\n",
    "I_n = np.eye(n)\n",
    "\n",
    "#create origin, destination, and origin-destination weight matrices\n",
    "Wo = np.kron(w, I_n)\n",
    "Wd = np.kron(I_n, w)\n",
    "Ww = np.kron(w, w)\n",
    "\n",
    "#create N x N identity matrix\n",
    "I_N = np.eye(N)\n",
    "\n",
    "#set origin rho, dest rho and o-d rho values\n",
    "o_rho = 0.5\n",
    "d_rho = 0.4\n",
    "od_rho = -0.2\n",
    "\n",
    "#compute jaobian term\n",
    "jac = I_N - (o_rho*Wo) - (d_rho*Wd) - (od_rho*Ww)\n",
    "\n",
    "#following the paper and code from LeSage to create matrices needd to compute effects\n",
    "imat_save = []\n",
    "dmat_save= []\n",
    "omat_save = []\n",
    "cnt = 0\n",
    "tt = np.arange(0,n)\n",
    "for j in tt:\n",
    "    imat = np.zeros((n,n))\n",
    "    dmat = np.zeros((n,n))\n",
    "    omat = np.zeros((n,n))\n",
    "    \n",
    "    dmat[cnt, :] = 1\n",
    "    dmat_save.append(dmat)\n",
    "    \n",
    "    omat[:, cnt] = 1\n",
    "    omat_save.append(omat)\n",
    "    \n",
    "    imat[cnt, cnt] = 1\n",
    "    imat_save.append(imat)\n",
    "    \n",
    "    cnt += 1\n",
    "\n",
    "#just converting the matrices to numpy array data types here with the required n^2 by n shape\n",
    "imat_save = np.array(imat_save).reshape((N,n))\n",
    "dmat_save = np.array(dmat_save).reshape((N,n))\n",
    "omat_save = np.array(omat_save).reshape((N,n))\n",
    "\n",
    "#set coefficient values\n",
    "beta_d = 1.\n",
    "beta_o = -0.5\n",
    "\n",
    "#create base matrices of effects\n",
    "tmp = np.zeros((N,n))\n",
    "tmp[dmat_save == 1] = beta_d\n",
    "tmp[omat_save == 1] = beta_o\n",
    "tmp[imat_save == 1] = beta_d + beta_o\n",
    " \n",
    "#compute total effects using jacobian term\n",
    "tmp = np.linalg.lstsq(jac, tmp)[0]\n",
    "total = (1./64.) * np.dot(np.dot(np.ones((1,64)), tmp), np.ones((8,1)))\n",
    "\n",
    "#compute intrazonal effects\n",
    "tmpi = np.zeros((N,n))\n",
    "tmpi[imat_save == 1] = tmp[imat_save == 1]\n",
    "intra = (1./64.) * np.dot(np.dot(np.ones((1,64)), tmpi), np.ones((8,1)))\n",
    "\n",
    "#compute origin effects\n",
    "tmpo = np.zeros((N,n))\n",
    "tmpo[omat_save == 1] = tmp[omat_save == 1]\n",
    "tmpo[imat_save == 1] = 0\n",
    "orig = (1./64.) * np.dot(np.dot(np.ones((1,64)), tmpo), np.ones((8,1)))\n",
    "\n",
    "#compute destination effects\n",
    "tmpd = np.zeros((N,n))\n",
    "tmpd[dmat_save == 1] = tmp[dmat_save == 1]\n",
    "tmpd[imat_save == 1] = 0\n",
    "dest = (1./64.) * np.dot(np.dot(np.ones((1,64)), tmpd), np.ones((8,1)))\n",
    "\n",
    "#compute network effects\n",
    "net = total - dest - orig - intra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "origin effects:  -0.621660531762\n",
      "dest effects:  1.21780957818\n",
      "intra effects:  0.0636121233569\n",
      "network effects:  1.05604928045\n",
      "total effects:  1.71581045023\n"
     ]
    }
   ],
   "source": [
    "#Lets check the compute effects here - presto! They match those from the paper (p. 200)\n",
    "print 'origin effects: ', orig[0][0]\n",
    "print 'dest effects: ', dest[0][0]\n",
    "print 'intra effects: ', intra[0][0]\n",
    "print 'network effects: ', net[0][0]\n",
    "print 'total effects: ', total[0][0]\n"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}