text: # -*- coding: utf-8 -*-
  """trying out VAR filtering and multidimensional fft
  
  Note: second half is copy and paste and doesn't run as script
  incomplete definitions of variables, some I created in shell
  
  Created on Thu Jan 07 12:23:40 2010
  
  Author: josef-pktd
  
  update 2010-10-22
  2 arrays were not defined, copied from fft_filter.log.py but I didn't check
  what the results are.
  Runs now without raising exception
  """
  
  import numpy as np
  from numpy.testing import assert_equal
  from scipy import signal
  from scipy.signal.signaltools import _centered as trim_centered
  _centered = trim_centered
  
  
  x = np.arange(40).reshape((2,20)).T
  x = np.arange(60).reshape((3,20)).T
  a3f = np.array([[[0.5,  1.], [1.,  0.5]],
                 [[0.5,  1.], [1.,  0.5]]])
  a3f = np.ones((2,3,3))
  
  
  nlags = a3f.shape[0]
  ntrim = nlags//2
  
  y0 = signal.convolve(x,a3f[:,:,0], mode='valid')
  y1 = signal.convolve(x,a3f[:,:,1], mode='valid')
  yf = signal.convolve(x[:,:,None],a3f)
  y = yf[:,1,:]  #
  yvalid = yf[ntrim:-ntrim,yf.shape[1]//2,:]
  #same result with fftconvolve
  #signal.fftconvolve(x[:,:,None],a3f).shape
  #signal.fftconvolve(x[:,:,None],a3f)[:,1,:]
  print trim_centered(y, x.shape)
  # this raises an exception:
  #print trim_centered(yf, (x.shape).shape)
  assert_equal(yvalid[:,0], y0.ravel())
  assert_equal(yvalid[:,1], y1.ravel())
  
  from scikits.statsmodels.tsa.filters import arfilter
  #copied/moved to statsmodels.tsa.filters
  def arfilter_old(x, a):
      '''apply an autoregressive filter to a series x
  
      x can be 2d, a can be 1d, 2d, or 3d
  
      Parameters
      ----------
      x : array_like
          data array, 1d or 2d, if 2d then observations in rows
      a : array_like
          autoregressive filter coefficients, ar lag polynomial
          see Notes
  
      Returns
      -------
      y : ndarray, 2d
          filtered array, number of columns determined by x and a
  
      Notes
      -----
  
      In general form this uses the linear filter ::
  
          y = a(L)x
  
      where
      x : nobs, nvars
      a : nlags, nvars, npoly
  
      Depending on the shape and dimension of a this uses different
      Lag polynomial arrays
  
      case 1 : a is 1d or (nlags,1)
          one lag polynomial is applied to all variables (columns of x)
      case 2 : a is 2d, (nlags, nvars)
          each series is independently filtered with its own
          lag polynomial, uses loop over nvar
      case 3 : a is 3d, (nlags, nvars, npoly)
          the ith column of the output array is given by the linear filter
          defined by the 2d array a[:,:,i], i.e. ::
  
              y[:,i] = a(.,.,i)(L) * x
              y[t,i] = sum_p sum_j a(p,j,i)*x(t-p,j)
                       for p = 0,...nlags-1, j = 0,...nvars-1,
                       for all t >= nlags
  
  
      Note: maybe convert to axis=1, Not
  
      TODO: initial conditions
  
      '''
      x = np.asarray(x)
      a = np.asarray(a)
      if x.ndim == 1:
          x = x[:,None]
      if x.ndim > 2:
          raise ValueError('x array has to be 1d or 2d')
      nvar = x.shape[1]
      nlags = a.shape[0]
      ntrim = nlags//2
      # for x is 2d with ncols >1
  
      if a.ndim == 1:
          # case: identical ar filter (lag polynomial)
          return signal.convolve(x, a[:,None], mode='valid')
          # alternative:
          #return signal.lfilter(a,[1],x.astype(float),axis=0)
      elif a.ndim == 2:
          if min(a.shape) == 1:
              # case: identical ar filter (lag polynomial)
              return signal.convolve(x, a, mode='valid')
  
          # case: independent ar
          #(a bit like recserar in gauss, but no x yet)
          result = np.zeros((x.shape[0]-nlags+1, nvar))
          for i in range(nvar):
              # could also use np.convolve, but easier for swiching to fft
              result[:,i] = signal.convolve(x[:,i], a[:,i], mode='valid')
          return result
  
      elif a.ndim == 3:
          # case: vector autoregressive with lag matrices
  #        #not necessary:
  #        if np.any(a.shape[1:] != nvar):
  #            raise ValueError('if 3d shape of a has to be (nobs,nvar,nvar)')
          yf = signal.convolve(x[:,:,None], a)
          yvalid = yf[ntrim:-ntrim, yf.shape[1]//2,:]
          return yvalid
  
  a3f = np.ones((2,3,3))
  y0ar = arfilter(x,a3f[:,:,0])
  print y0ar, x[1:] + x[:-1]
  yres = arfilter(x,a3f[:,:,:2])
  print np.all(yres == (x[1:,:].sum(1) + x[:-1].sum(1))[:,None])
  
  # don't do these imports, here just for copied fftconvolve
  from scipy.fftpack import fft, ifft, ifftshift, fft2, ifft2, fftn, \
       ifftn, fftfreq
  from numpy import product,array
  
  from scikits.statsmodels.tsa.filters.filtertools import fftconvolveinv as fftconvolve
  #copied/moved to statsmodels.tsa.filters
  def fftconvolve_old(in1, in2, in3=None, mode="full"):
      """Convolve two N-dimensional arrays using FFT. See convolve.
  
      copied from scipy.signal.signaltools, but here used to try out inverse filter
      doesn't work or I can't get it to work
  
       2010-10-23:
      looks ok to me for 1d,
      from results below with padded data array (fftp)
      but it doesn't work for multidimensional inverse filter (fftn)
      original signal.fftconvolve also uses fftn
  
      """
      s1 = array(in1.shape)
      s2 = array(in2.shape)
      complex_result = (np.issubdtype(in1.dtype, np.complex) or
                        np.issubdtype(in2.dtype, np.complex))
      size = s1+s2-1
  
      # Always use 2**n-sized FFT
      fsize = 2**np.ceil(np.log2(size))
      IN1 = fftn(in1,fsize)
      #IN1 *= fftn(in2,fsize) #JP: this looks like the only change I made
      IN1 /= fftn(in2,fsize)  # use inverse filter
      # note the inverse is elementwise not matrix inverse
      # is this correct, NO  doesn't seem to work for VARMA
      fslice = tuple([slice(0, int(sz)) for sz in size])
      ret = ifftn(IN1)[fslice].copy()
      del IN1
      if not complex_result:
          ret = ret.real
      if mode == "full":
          return ret
      elif mode == "same":
          if product(s1,axis=0) > product(s2,axis=0):
              osize = s1
          else:
              osize = s2
          return _centered(ret,osize)
      elif mode == "valid":
          return _centered(ret,abs(s2-s1)+1)
  
  yff = fftconvolve(x.astype(float)[:,:,None],a3f)
  
  rvs = np.random.randn(500)
  ar1fft = fftconvolve(rvs,np.array([1,-0.8]))
  #ar1fftp = fftconvolve(np.r_[np.zeros(100),rvs,np.zeros(100)],np.array([1,-0.8]))
  ar1fftp = fftconvolve(np.r_[np.zeros(100),rvs],np.array([1,-0.8]))
  ar1lf = signal.lfilter([1], [1,-0.8], rvs)
  
  ar1 = np.zeros(501)
  for i in range(1,501):
      ar1[i] = 0.8*ar1[i-1] + rvs[i-1]
  
  #the previous looks wrong, is for generating ar with delayed error,
  #or maybe for an ma(1) filter, (generating ar and applying ma filter are the same)
  #maybe not since it replicates lfilter and fftp
  #still strange explanation for convolution
  #ok. because this is my fftconvolve, which is an inverse filter (read the namespace!)
  
  #This is an AR filter
  errar1 = np.zeros(501)
  for i in range(1,500):
      errar1[i] = rvs[i] - 0.8*rvs[i-1]
  
  #print ar1[-10:]
  #print ar1fft[-11:-1]
  #print ar1lf[-10:]
  #print ar1[:10]
  #print ar1fft[1:11]
  #print ar1lf[:10]
  #print ar1[100:110]
  #print ar1fft[100:110]
  #print ar1lf[100:110]
  #
  #arloop - lfilter - fftp (padded)  are the same
  print '\n compare: \nerrloop - arloop - fft - lfilter - fftp (padded)'
  #print np.column_stack((ar1[1:31],ar1fft[:30], ar1lf[:30]))
  print np.column_stack((errar1[1:31], ar1[1:31],ar1fft[:30], ar1lf[:30],
                         ar1fftp[100:130]))
  
  def maxabs(x,y):
      return np.max(np.abs(x-y))
  
  print maxabs(ar1[1:], ar1lf)  #0
  print maxabs(ar1[1:], ar1fftp[100:-1]) # around 1e-15
  
  rvs3 = np.random.randn(500,3)
  a3n = np.array([[1,1,1],[-0.8,0.5,0.1]])
  a3n = np.array([[1,1,1],[-0.8,0.0,0.0]])
  a3n = np.array([[1,-1,-1],[-0.8,0.0,0.0]])
  a3n = np.array([[1,0,0],[-0.8,0.0,0.0]])
  a3ne = np.r_[np.ones((1,3)),-0.8*np.eye(3)]
  a3ne = np.r_[np.ones((1,3)),-0.8*np.eye(3)]
  ar13fft = fftconvolve(rvs3,a3n)
  
  ar13 = np.zeros((501,3))
  for i in range(1,501):
      ar13[i] = np.sum(a3n[1,:]*ar13[i-1]) + rvs[i-1]
  
  #changes imp was not defined, not sure what it is supposed to be
  #copied from a .log file
  imp = np.zeros((10,3))
  imp[0]=1
  
  a3n = np.array([[1,0,0],[-0.8,0.0,0.0]])
  fftconvolve(np.r_[np.zeros((100,3)),imp],a3n)[100:]
  a3n = np.array([[1,0,0],[-0.8,-0.50,0.0]])
  fftconvolve(np.r_[np.zeros((100,3)),imp],a3n)[100:]
  
  a3n3 = np.array([[[ 1. ,  0. ,  0. ],
                   [ 0. ,  1. ,  0. ],
                   [ 0. ,  0. ,  1. ]],
  
                  [[-0.8,  0. ,  0. ],
                   [ 0. , -0.8,  0. ],
                   [ 0. ,  0. , -0.8]]])
  
  a3n3 = np.array([[[ 1. ,  0.5 ,  0. ],
                    [ 0. ,  1. ,  0. ],
                    [ 0. ,  0. ,  1. ]],
  
                   [[-0.8,  0. ,  0. ],
                    [ 0. , -0.8,  0. ],
                    [ 0. ,  0. , -0.8]]])
  ttt = fftconvolve(np.r_[np.zeros((100,3)),imp][:,:,None],a3n3.T)[100:]
  gftt = ttt/ttt[0,:,:]
  
  a3n3 = np.array([[[ 1. ,  0 ,  0. ],
                    [ 0. ,  1. ,  0. ],
                    [ 0. ,  0. ,  1. ]],
  
                   [[-0.8,  0.2 ,  0. ],
                    [ 0 ,  0.0,  0. ],
                    [ 0. ,  0. , 0.8]]])
  ttt = fftconvolve(np.r_[np.zeros((100,3)),imp][:,:,None],a3n3)[100:]
  gftt = ttt/ttt[0,:,:]
  signal.fftconvolve(np.dstack((imp,imp,imp)),a3n3)[1,:,:]
  
  nobs = 10
  imp = np.zeros((nobs,3))
  imp[1] = 1.
  ar13 = np.zeros((nobs+1,3))
  for i in range(1,nobs+1):
      ar13[i] = np.dot(a3n3[1,:,:],ar13[i-1]) + imp[i-1]
  
  a3n3inv = np.zeros((nobs+1,3,3))
  a3n3inv[0,:,:] = a3n3[0]
  a3n3inv[1,:,:] = -a3n3[1]
  for i in range(2,nobs+1):
      a3n3inv[i,:,:] = np.dot(-a3n3[1],a3n3inv[i-1,:,:])
  
  
  a3n3sy = np.array([[[ 1. ,  0 ,  0. ],
                    [ 0. ,  1. ,  0. ],
                    [ 0. ,  0. ,  1. ]],
  
                   [[-0.8,  0.2 ,  0. ],
                    [ 0 ,  0.0,  0. ],
                    [ 0. ,  0. , 0.8]]])
  
  nobs = 10
  a = np.array([[[ 1. ,  0. ],
          [ 0. ,  1. ]],
  
         [[-0.8,  0.0 ],
          [ -0.1 , -0.8]]])
  
  
  a2n3inv = np.zeros((nobs+1,2,2))
  a2n3inv[0,:,:] = a[0]
  a2n3inv[1,:,:] = -a[1]
  for i in range(2,nobs+1):
      a2n3inv[i,:,:] = np.dot(-a[1],a2n3inv[i-1,:,:])
  
  nobs = 10
  imp = np.zeros((nobs,2))
  imp[0,0] = 1.
  
  #a2 was missing, copied from .log file, not sure if correct
  a2 = np.array([[[ 1. ,  0. ],
          [ 0. ,  1. ]],
  
         [[-0.8,  0. ],
          [0.1, -0.8]]])
  
  ar12 = np.zeros((nobs+1,2))
  for i in range(1,nobs+1):
      ar12[i] = np.dot(-a2[1,:,:],ar12[i-1]) + imp[i-1]
  
  u = np.random.randn(10,2)
  ar12r = np.zeros((nobs+1,2))
  for i in range(1,nobs+1):
      ar12r[i] = np.dot(-a2[1,:,:],ar12r[i-1]) + u[i-1]
  
  a2inv = np.zeros((nobs+1,2,2))
  a2inv[0,:,:] = a2[0]
  a2inv[1,:,:] = -a2[1]
  for i in range(2,nobs+1):
      a2inv[i,:,:] = np.dot(-a2[1],a2inv[i-1,:,:])
  
  import scipy.stats as stats
  import numpy as np
  nbins = 12
  binProb = np.zeros(nbins) + 1.0/nbins
  binSumProb = np.add.accumulate(binProb)
  print binSumProb
  print stats.gamma.ppf(binSumProb,0.6379,loc=1.6,scale=39.555)
  