Permalink
Browse files

sandbox

tsa/arima small bugfix in currently unused function and add arma_pacf
add normal jacobian to regression tools
  • Loading branch information...
1 parent de54b07 commit 6fa6b94c9f80bbe58696ee95a2a4976cac730747 @josef-pkt josef-pkt committed Mar 11, 2010
Showing with 122 additions and 1 deletion.
  1. +2 −0 .bzrignore
  2. +101 −0 scikits/statsmodels/sandbox/regression/tools.py
  3. +19 −1 scikits/statsmodels/sandbox/tsa/arima.py
View
@@ -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
@@ -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)
+
+
+
+
@@ -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

0 comments on commit 6fa6b94

Please sign in to comment.