sandbox

tsa/arima small bugfix in currently unused function and add arma_pacf
add normal jacobian to regression tools

.bzrignore

@@ -19,9 +19,11 @@ build
 *.diff
 .settings/
 *.svn/
+*.log.py
 # Egg metadata
 ./*.egg-info
 # The shelf plugin uses this dir
 ./.shelf
 # Mac droppings
 .DS_Store
+help

scikits/statsmodels/sandbox/regression/tools.py

@@ -0,0 +1,101 @@
+'''gradient/Jacobian of normal loglikelihood
+
+use chain rule
+
+
+A: josef-pktd
+
+'''
+
+
+import numpy as np
+
+
+def norm_lls(y, params):
+    '''normal loglikelihood given observations and mean mu and variance sigma2
+
+    Parameters
+    ----------
+    y : array, 1d
+        normally distributed random variable
+    params: array, (nobs, 2)
+        array of mean, variance (mu, sigma2) with observations in rows
+    '''
+
+    mu, sigma2 = params.T
+    lls = -0.5*(np.log(2*np.pi) + np.log(sigma2) + (y-mu)**2/sigma2)
+    return lls
+
+def norm_lls_grad(y, params):
+    '''Jacobian of normal loglikelihood wrt mean mu and variance sigma2
+
+    Parameters
+    ----------
+    y : array, 1d
+        normally distributed random variable
+    params: array, (nobs, 2)
+        array of mean, variance (mu, sigma2) with observations in rows
+
+    Returns
+    -------
+    grad : array (nobs, 2)
+        derivative of loglikelihood for each observation wrt mean in first
+        column, and wrt variance in second column
+    '''
+    mu, sigma2 = params.T
+    dllsdmu = (y-mu)/sigma2
+    dllsdsigma2 = ((y-mu)**2/sigma2 - 1)/np.sqrt(sigma2)
+    return np.column_stack((dllsdmu, dllsdsigma2))
+
+
+def mean_grad(x, beta):
+    '''gradient/Jacobian for d (x*beta)/ d beta
+    '''
+    return x
+
+def normgrad(y, x, params):
+    '''Jacobian of normal loglikelihood wrt mean mu and variance sigma2
+
+    Parameters
+    ----------
+    y : array, 1d
+        normally distributed random variable with mean x*beta, and variance sigma2
+    x : array, 2d
+        explanatory variables, observation in rows, variables in columns
+    params: array_like, (nvars + 1)
+        array of coefficients and variance (beta, sigma2)
+
+    Returns
+    -------
+    grad : array (nobs, 2)
+        derivative of loglikelihood for each observation wrt mean in first
+        column, and wrt variance in second column
+    assume params = (beta, sigma2)
+
+    Notes
+    -----
+    TODO: for heteroscedasticity need sigma to be a 1d array
+
+    '''
+    beta = params[:-1]
+    sigma2 = params[-1]*np.ones((len(y),1))
+    dmudbeta = mean_grad(x, beta)
+    mu = np.dot(x, beta)
+    #print beta, sigma2
+    params2 = np.column_stack((mu,sigma2)) #np.concatenate((beta, np.atleast_1d(sigma2)))
+    dllsdms = norm_lls_grad(y,params2)
+    grad = np.column_stack((dllsdms[:,:1]*dmudbeta, dllsdms[:,:1]))
+    return grad
+
+sig = 0.1
+beta = np.ones(2)
+rvs = np.random.randn(10,3)
+x = rvs[:,1:]
+y = np.dot(x,beta) + sig*rvs[:,0]
+
+params = [1,1,0.1**2]
+print normgrad(y, x, params)
+
+
+
+

scikits/statsmodels/sandbox/tsa/arima.py

@@ -257,6 +257,24 @@ def arma_acf(ar, ma, nobs=10):
     acovf = arma_acovf(ar, ma, nobs)
     return acovf/acovf[0]
 
+def arma_pacf(ar, ma, nobs=10):
+    '''partial autocorrelation function of an ARMA process
+
+    Notes
+    -----
+    solves yule-walker equation for each lag order up to nobs lags
+
+    not tested/checked yet
+    '''
+    apacf = np.zeros(nobs)
+    acov = tsa.arma_acf(ar,ma)
+
+    apacf[0] = 1.
+    for k in range(2,nobs+1):
+        r = acov[:k];
+        apacf[k] = np.linalg.solve(scipy.linalg.toeplitz(r[:-1]), r[1:])[-1]
+
+
 def arma_impulse_response(ar, ma, nobs=100):
     '''get the impulse response function (MA representation) for ARMA process
 
@@ -365,7 +383,7 @@ def lpol2index(ar):
         index (lags) of lagpolynomial with non-zero elements
     '''
 
-    index = np.nonzero(ar)
+    index = np.nonzero(ar)[0]
     coeffs = ar[index]
     return coeffs, index