pyprocess.py

from __future__ import division

"""
PyProcess
@author: Cameron Davidson-Pilon
"""
import numpy as np
import scipy as sp
import scipy.stats as stats
from scipy.special import gammainc
import warnings
import matplotlib.pyplot as plt
import pdb

np.seterr( divide="raise")


#TODO change all generate_ to sample_



class Step_process(object):
    """
      This is the class of finite activity jump/event processes (eg poisson process, renewel process etc). 
      Parameters:
        time_space_constraints: a dictionary contraining AT LEAST "startTime":float, "startPosition":float,  

    """

    def  __init__(self, startTime = 0, startPosition = 0, endTime = None, endPosition = None):
    
        self.conditional=False
        self.startTime = startTime
        self.startPosition = startPosition
        if ( endTime != None ) and ( endPosition != None ):
            assert endTime > startTime, "invalid parameter: endTime > startTime"
            self.endTime = endTime
            self.endPosition = endPosition
            self.condition = True
        elif ( endTime != None ) != ( endPosition != None ):
            raise Exception( "invalid parameter:", "Must include both endTime AND endPosition or neither" )

            
    def _check_time(self,t):
        if t<self.startTime:
            raise Exception("Attn: inputed time not valid. Check if beginning time is less than startTime (0 by default).")
            
                
    def sample_path(self,times, N=1):
        try:
            self._check_time(times[0])
        except TypeError:
            self._check_time( times )
            
        return self._sample_path(times, N)
    
    
    def sample_jumps(self,T, N = 1):
        self._check_time(T)
        return self._sample_jumps(T,N)
    
    def mean(self,t):
        self._check_time(t)
        return self._mean(t)
    
    def var(self,t):
        self._check_time(t)
        return self._var(t)
    
    def mean_number_of_jumps(self,t):
        self._check_time(t)
        return self._mean_number_of_jumps(t)
    
    def sample_position(self,t, N=1):
        self._check_time(t)
        return self._sample_position(t, N)
    
    
    def _mean_number_of_jumps(self,t, n_simulations = 10000):
        self._check_time(t)
        """
        This needs to solve the following:
          E[ min(N) | T_1 + T_2 + .. T_N >= t ]
        <=>
          E[ N | T_1 + T_2 + .. T_{N-1} < t, T_1 + T_2 + .. T_N >= t ]
        where we can sample T only. 
        
        """
        warnings.warn("Attn: performing MC simulation. %d simulations."%n_simulations)
        v = np.zeros( n_simulations )
        continue_ = True
        i = 1
        t = t - self.startTime #should be positive.
        size = v.shape[0]
        sum = 0
        sum_sq = 0
        while continue_:
            v += self.T.rvs( size )
            #count how many are greater than t
            n_greater_than_t = (v>=t).sum()
            sum += n_greater_than_t*i 
            
            v = v[v < t]
            i+=1
            size = v.shape[0]
            empty = size == 0
            
            
            if empty:
                continue_ = False
                
            if i > 10000:
                warnings.warn("The algorithm did not converge. Be sure 'T' is a non negative random \
                                variable, or set 't' lower.")
                break

        return sum/n_simulations
        
    def _var(self,t):
        self._check_time(t)
        return self._mean_number_of_jumps(t)*self.J.var()


    def plot(self, t,N=1, **kwargs):
        assert N >= 1, "N must be greater than 0."
        for n in range(N):
            path = [self.startPosition, self.startPosition]
            times = [self.startTime]
            _path, _times = self.sample_path(t,1)
            
            for i in range(2*_path.shape[1]):
                j = int( (i - i%2)/2)
                path.append( _path[0][j] )
                times.append( _times[0][j] )
            #path.append( _path[0][-1] )
            times.append( t) 
            #pdb.set_trace()
            plt.plot( times, path, **kwargs )
        plt.xlabel( "time, $t$")
        plt.ylabel( "position of process")
        plt.show()
        return
        

        
class Renewal_process(Step_process):
    """    
    parameters:
        T: T is a scipy "frozen" random variable object, from the scipy.stats library, that has the same distribution as the inter-arrival times 
    i.e. t_{i+1} - t_{i} equal in distribution to T, where t_i are jump/event times.
    T must be a non-negative random variable, possibly constant. The constant case in included in this library, it 
    can be called as Contant(c), where c is the interarrival time.
    
        J: is a scipy "frozen" random variable object that has the same distribution as the jump distribution. It can have any real support. 
    Ex: for a poisson process, J is equal to 1 with probability 1, and we can use the Constant(1) class, 
    but J is random for the compound poisson process.
    
        Note 1. endTime and endPosition are not allowed for this class.
        Note 2. If you are interested in Poisson or Compound poisson process, they are implemented in PyProcess.
                Those implementations are recommened to use versus using a Renewal_process. See Poisson_process
                and Compound_poisson_process.
    Ex:
    import scipy.stats as stats
    
    #create a standard renewal process
    renewal_ps= Renewal_process( J = Constant(3.5), T = stats.poisson(1), startTime = 0, startPosition = 0 )
    
    """
    
    
    def __init__(self, T, J, startPosition = 0, startTime = 0):
    
        super(Renewal_process, self).__init__(startTime = startTime,startPosition = startPosition)
        self.T = T
        self.J = J
        self.renewal_rate = 1/self.T.mean()

    def forward_recurrence_time_pdf(self,x):
        "the forward recurrence time RV is the time need to wait till the next event, after arriving to the system at a large future time"
        return self.renewal_rate*self.T.sf(x)
    
    def backwards_recurrence_time_pdf(self,x):
        "the backwards recurrence time is the time since the last event after arriving to the system at a large future time"
        return self.forward_recurrence_time_pdf(x)
    
    def forward_recurrence_time_cdf(self,x):
        "the forward recurrence time RV is the time need to wait till the next event, after arriving to the system at a large future time"
        return int.quad(self.forward_recurrence_time_pdf, 0, x)
    
    def backward_recurrence_time_cdf(self,x):
        "the backwards recurrence time is the time since the last event after arriving to the system at a large future time"
        return int.quad(self.forward_recurrence_time_pdf, 0, x) 
    
    def spread_pdf(self,x):
        """the RV distributed according to the spread_pdf is the (random) length of time between the previous and next event/jump
         when you arrive at a large future time."""
        try:
            return self.renewal_rate*x*self.T.pdf(x)
        except AttributeError:
            return self.renewal_rate*x*self.T.pmf(x)
    
    def _sample_path(self,time, N =1):
        """
        parameters: 
            time:  an interable that returns floats OR if a single float, returns (path,jumps) up to time t.
            N: the number of paths to return. negative N returns an iterator of size N.
        returns:
            path: the sample path at the times in times
            jumps: the times when the process jumps
        
        Scales linearly with number of sequences and time.
    
        """
        if N < 0:
            return self._iterator_sample_paths(time, -N)
        
        else:
            return self.__sample_path( time, N ) #np.asarray( [self.__sample_path(times) for i in xrange(N) ] )
        
    def _iterator_sample_paths( self, times, N):
        
        for i in xrange(N):
            yield self.__sample_path(times) 
    
    def __sample_path( self, times, N=1 ):
        
        if isinstance( times, int): #not the best way to do this.
        
            """user wants a path up to time times."""
            def generate_jumps( x, J, start_position):
                n = x.shape[0]
                return start_position + J.rvs(n).cumsum()
                
            _times = self.sample_jumps(times,N)
            _path =  np.asarray( map( lambda x: generate_jumps(x, self.J, self.startPosition), _times) ) 
            return (_path, _times )
            
        else:
           # times is an iterable, user requests "give me positions at these times."
           pass
           #TODO
        
        
       
        
    
    def __sample_jumps(self,  t ):
        
        quantity_per_i = 10
        times = np.array([])
        c = 0 #measure of size
        while True:
            c += quantity_per_i
            times = np.append( times, self.T.rvs(quantity_per_i) )
            times = times[ times.cumsum() < t ]
            if times.shape[0] < c:
                #this will only happen if the list is truncated. 
                return times.cumsum()

    def iterator_sample_jumps(self, t, N):
    
        for i in xrange(N):
            yield __sample_jumps(t)

        
    def _sample_jumps(self,t, N=1):
        """        
        Generate N simulations of jump times prior to some time t > startTime.
        
        parameters:
            t: the end of the processes. 
            N: the number of samples to return, -N to return a generator of size N
        returns:
            An irregular numpy array, with first dimension N, and each element an np.array 
            with random length.
            
        example:
        >> print renewal_ps.generate_sample_jumps( 1, N=2)
        np.array( [ [2,2.5,3.0], [1.0,2.0] ], dtype=0bject )
        
        """
        if N<0:
            return (  self.__sample_jumps(t) for i in xrange(-N) )
        else:
            return np.asarray( [ self.__sample_jumps(t) for i in xrange(N) ] )
    
    
    def _mean(self,t):
        #try:
            #return self.J.mean()*self.T.mean()   
        #except:
            return self.J.mean()*self.mean_number_of_jumps(t)     
    
    def _var(self,t):
            return self.mean_number_of_jumps(t)*self.J.var()
    
    def _sample_position(self,t, N = 1):
        """
        Generate the position of the process at time t > startTime
        parameters:
            N: number of samples
            t: position to sample at.
        returns:
            returns a (n,) np.array 
        """
        v = np.zeros( N ) 
        iv = np.ones( N, dtype=bool )
        x = self.startPosition*np.ones( N) 
        continue_ = True
        i = 1
        t = t - self.startTime #should be positive.
        size = v.shape[0]
        while continue_:
            n_samples = iv.sum()
            v[iv] += self.T.rvs( n_samples )
            
            #if the next time is beyond reach, we do not add a new jump
            iv[ v >= t] = False
            n_samples = iv.sum()

            x[iv] += self.J.rvs( n_samples )

            size = iv.sum() #any left?
            empty = size == 0
            
            i+=1
            if empty:
                continue_ = False
                
            if i > 10000:
                warnings.warn("The algorithm did not converge. Be sure 'T' is a non negative random \
                                variable, or set 't' to a smaller value.")
                break
                
        return x
            
class Poisson_process(Renewal_process):
    """
    This implements a Poisson process with rate parameter 'rate' (lambda was already taken!). 
    Parameters:
        rate: float > 0 that define the inter-arrival rate. 
        startTime: (default 0) The start time of the process
        startPosition: (default 0) The starting position of the process at startTime
        
        conditional processes:
          If the process is known at some future time, endPosition and endTime 
          condition the process to arrive there.
        endTime: (default None) the time in the future the process is known at, > startTime
        endPosition: (default None) the position the process is at in the future.
                     > startPosition and equal to startPosition + int.
         
    ex:
    
    #condition the process to be at position 5 at time 10.
    pp = Poisson_process( rate = 3, endTime = 10, endPosition = 5 )  
    
    """
    
    def __init__(self,rate=1, startTime = 0, startPosition = 0, endPosition = None, endTime = None):
        
       assert rate > 0, "invalid parameter: rate parameter must be greater than 0."
       self.rate = rate
       self.Exp = stats.expon(1/self.rate)
       self.Con = Constant(1) 
       super(Poisson_process,self).__init__(J=self.Con, T=self.Exp, startTime = startTime, startPosition = startPosition)
       self.Poi = stats.poisson
       if ( endTime != None ) and ( endPosition != None ):
            assert endTime > startTime, "invalid times: endTime > startTime."
            self.endTime = endTime
            self.endPosition = endPosition
            self.condition = True
            self.Bin = stats.binom
       elif  ( endTime != None ) != ( endPosition != None ):
            raise Exception( "invalid parameter:", "Must include both endTime AND endPosition or neither" )
       
    def _mean(self,t):
        """
        recall that a conditional poisson process N_t | N_T=n ~ Bin(n, t/T)
        """
        if not self.conditional:  
            return self.startPosition + self.rate*(t-self.startTime)
        else:
            return self.endPosition*float(t)/self.endTime
    
    def _var(self,t):
        """
        parameters:
            t: a time > startTime (and less than endTime if present).
        returns:
            a float.
            
        recall that a conditional poisson process N_t | N_T=n ~ Bin(n, t/T)
        
        """
        if self.conditional:
            return self.endPosition*(1-float(t)/self.endTime)*float(t)/self.endTime
        else:
            return self.rate*(t-self.startTime)
    
    def __sample_jumps(self,t):
        
        if self.conditional:
            p = self.Bin.rvs(self.endPosition, t/self.endTime)
        else:
            p = self.Poi.rvs(self.rate*(t-self.startTime))
        path=[]
        #array = [self.startTime+(T-self.startTime)*np.random.random() for i in xrange(p)]
        jump_times = self.startTime + (t - self.startTime)*np.random.random(p)
        jump_times.sort() 
        #for i in xrange(p):
        #    x+=1
        #    path.append((array[i],x))
        #    i+=1
        return  jump_times 
        
    def _sample_jumps(self, t, N=1):
        """
        Probably to be deleted.
            T:  a float
            N: the number of sample paths
        Returns:
            an (N,) np.array of jump times. 
        """
        if N<0:
            return ( self.__sample_jumps(t) for i in xrange(-N) )
        else:
            return np.asarray( [self.__sample_jumps(t) for i in xrange(N)] )
    
    
    
    def _sample_position(self,t, N=1):
        if self.conditional:
            return self.Bin.rvs(self.endPosition-self.startPosition, float(t)/self.endTime, size = N)+self.startPosition
        else:
            return self.Poi.rvs(self.rate*(t-self.startTime), size = N)+self.startPosition

            
class Marked_poisson_process(Renewal_process):
    """
    This class constructs marked poisson process i.e. at exponentially distributed times, a 
    Uniform(L,U) is generated.
    
    parameters:
        rate: the rate for the poisson process
        U: the upper bound of uniform random variate generation
        L: the lower bound of uniform random variate generation
        
    
    Note there are no other time-space constraints besides startTime
    """
    def __init__(self,rate =1, L= 0, U = 1, startTime = 0):
        self.Poi = stats.poisson
        self.L = L
        self.U = U
        self.startTime = startTime
        self.rate = rate

        
    def generate_marked_process(self,T):
                
        p = self.Poi.rvs(self.rate*(T-self.startTime))

        times = self.startTime + (T-self.startTime)*np.random.random( p )
        path = self.L+(self.U-self.L)*np.random.random( p )
        
        return ( path, times )
        

class Compound_poisson_process(Renewal_process):
    """
    This process has expontially distributed inter-arrival times (i.e. 'rate'-poisson distributed number of jumps at any time),
    and has jump distribution J.
    parameters:
    
     J: a frozen scipy.stats random variable instance. It can have any support.
     
     Note: endTime and endPosition constraints will not be statisfied.
     
     Example:
     >> import stats.scipy as stats
     
     >> Nor = stats.norm(0,1)
     >> cmp = Compound_poisson_process(J = Nor) 
    
    """
    def __init__(self, J, rate = 1,startTime = 0, startPosition = 0):      
        assert rate > 0, "Choose rate to be greater than 0."
        self.J = J
        self.rate = rate
        self.Exp = stats.expon(1./self.rate)
        super(Compound_poisson_process, self).__init__(J = J, T= T, startTime = startTime, startPosition = startPosition)
        self.Poi = stats.poisson
        
    def _mean(self,t):
        return self.startPosition + self.rate*(t-self.startTime)*self.J.mean()
    
    def _var(self,t):
        return self.rate*(t-self.startTime)*(self.J.var()-self.J.mean()**2)
    """ 
    def __sample_jumps(self,T):
        p = self.Poi.rvs(self.rate*(T-self.startTime))
        x = self.startPosition
        path=[]
        array = [self.startTime+(T-self.startTime)*np.random.rand() for i in range(p)]
        array.sort() 
        for i in range(p):
            x+=self.J.rvs()
            path.append((array[i],x))
            i+=1
        return path
    """


class Diffusion_process(object):
    #
    #    Class that can be overwritten in the subclasses:
    #        _sample_position(t, N=1)
    #        _mean(t)
    #        _var(t)
    #        _sample_path(times, N=1)
    #
    #    Class that should be present in subclasses:
    #
    #        _transition_pdf(x,t,y)
    #
    #
    
    def  __init__(self, startTime = 0, startPosition = 0, endTime = None, endPosition = None):
    
        self.conditional=False
        self.startTime = startTime
        self.startPosition = startPosition
        if ( endTime != None ) and ( endPosition != None ):
            assert endTime > startTime, "invalid parameter: endTime > startTime"
            self.endTime = endTime
            self.endPosition = endPosition
            self.conditional = True
        elif ( endTime != None ) != ( endPosition != None ):
            raise Exception( "invalid parameter:", "Must include both endTime AND endPosition or neither" )

    
    def transition_pdf(self,t,y):
        self._check_time(t)
        "this method calls self._transition_pdf(x,t,y) in the subclass"
        try:
            if not self.conditional:
                return self._transition_pdf(self.startPosition, t-self.startTime, y)
            else:
                return self._transition_pdf(self.startPosition, t-self.startTime, y)*self._transition_pdf(y, self.endTime-t, self.endPosition)\
                        /self._transition_pdf(self.startPosition,self.endTime - self.startTime, self.endPosition)
        except AttributeError:
            raise AttributeError("Attn: transition density for process is not defined.")

    def expected_value(self,t, f= lambda x:x, N=1e6):
        """
        This function calculates the expected value of  E[ X_t | F ] where F includes start conditions and possibly end conditions.
        
        
        """
        warnings.warn( "Performing Monte Carlo with %d simulations."%N)
        self._check_time(t)
        if not self.conditional:
            return f( self.sample_position(t, N) ).mean() 
        else:
            #This uses a change of measure technique.
            """
            sum=0
            self.conditional=False
            for i in range(N):
                X = self.generate_position_at(t)
                sum+=self._transition_pdf(X,self.endTime-t,self.endPosition)*f(X)
            self.conditional=True
            """
            x = self.sample_position(t, N)
            mean = (self._transition_pdf(x,self.endTime-t,self.endPosition)*f(x)).mean()
            return mean/self._transition_pdf(self.startPosition, self.endTime-self.startTime, self.endPosition)
        
    def sample_position(self,t, N=1):
        """        
        if _get_position_at() is not overwritten in a subclass, this function will use euler scheme
        """

        self._check_time(t)
        return self._sample_position(t, N)
        
    
    def mean(self,t):
        self._check_time(t)
        return self._mean(t)
    
    
    def var(self,t):
        self._check_time(t)
        return self._var(t)
    
    def sample_path(self,times, N = 1):
        self._check_time(times[0])
        return self._sample_path(times, N)
    
    
    
    
    
    def _sample_path(self,times, N=1 ):
        return self.Euler_scheme(times, N)
    
    def _var(self,t, N = 1e6):
        """
        var = SampleVarStat()
        for i in range(10000):
            var.push(self.generate_position_at(t))
        return var.get_variance()
        """
        return self.sample_position(t, n).var()
    
    def _mean(self,t):
        return self.expected_value( t)
    
    def _sample_position(self,t):
        return self.Euler_scheme(t)
        
    def _transition_pdf(self,x,t,y):
        warning.warn( "Attn: transition pdf not defined" )
    
    def _check_time(self,t):
        if t<self.startTime:
            warnings.warn( "Attn: inputed time not valid (check if beginning time is less than startTime)." )
    
    def Euler_scheme(self, times,delta=0.001):
        """
        times is an array of floats.
        The process needs the methods drift() and diffusion() defined.
        """
        warnings.warn("Attn: starting an Euler scheme to approxmiate process.")
        
        Nor = stats.norm()
        finalTime = times[-1]
        steps = int(finalTime/delta)
        t = self.startTime
        x=self.startPosition
        path=[]
        j=0
        time = times[j]
        for i in xrange(steps):
            if t+delta>time>t:
                delta = time-t
                x += drift(x,t)*delta + np.sqrt(delta)*diffusion(x,t)*Nor.rvs()
                path.append((x,time))
                delta=0.001
                j+=1
                time = times[j]
            else:
                x += drift(x,t)*delta + np.sqrt(delta)*diffusion(x,t)*Nor.rvs()
                t += delta
            
        return path  

    def Milstein_Scheme(self, times, delta = 0.01 ):
        if not all( map( lambda x: hasattr(self, x), ( 'drift', 'diffusion', 'diffusion_prime' ) ) ):
            raise AttributeError("The process does not have 'drift', 'diffusion', or 'diffusion_prime' methods")
        
        pass
    
    
        
        
    def plot(self, times ,N=1, **kwargs):
        assert N >= 1, "N must be greater than 0."
        try:
            self._check_time(times[-1] )
            plt.plot(times, self.sample_path(times, N).T, **kwargs )
        except:
            self._check_time(times)
            times = np.linspace(self.startTime, times, 100)
            path = self.sample_path(times, N).T
            plt.plot(times, path, **kwargs )            
        plt.xlabel( "time, $t$")
        plt.ylabel( "position of process")
        plt.show()
        return
        
class Wiener_process(Diffusion_process):
    """
    This implements the famous Wiener process. I choose not to call it Brownian motion, as 
    brownian motion is a special case of this with 0 drift and variance equal to t.
    
    dW_t = mu*dt + sigma*dB_t
    
    W_t ~ N(mu*t, sigma**2t)
    
    parameters: 
        mu: the constant drift term, float
        sigma: the constant volatility term, float > 0
    """
    
    def __init__(self, mu, sigma, startTime = 0, startPosition = 0, endPosition = None, endTime = None):
        super(Wiener_process,self).__init__(startTime, startPosition, endTime, endPosition)
        self.mu = mu
        self.sigma = sigma
        self.Nor = stats.norm()
    
    def _transition_pdf(self,x,t,y):
        return np.exp(-(y-x-self.mu*(t-self.startTime))**2/(2*self.sigma**2*(t-self.startTime)))\
            /np.sqrt(2*pi*self.sigma*(t-self.startTime))    
    
    def _mean(self,t):
        if self.conditional:
            delta1 = t - self.startTime
            delta2 = self.endTime - self.startTime
            return self.startPosition + self.mu*delta1 + (self.endPosition-self.startPosition-self.mu*delta2)*delta1/delta2
        else:
            return self.startPosition+self.mu*(t-self.startTime)

    def _var(self,t):
        if self.conditional:
            delta1 = self.sigma**2*(t-self.startTime)*(self.endTime-t)
            delta2 = self.endTime-self.startTime
            return delta1/delta2
        else:
            return self.sigma**2*(t-self.startTime)
        
    def _sample_position(self,t, n=1):
        """
        This incorporates both conditional and unconditional
        """
        return self.mean(t) + np.sqrt(self.var(t))*self.Nor.rvs(n)

    def _sample_path(self,times, N = 1):

        path=np.zeros( (N,len(times)) )
        path[ :, 0] = self.startPosition
        times = np.insert( times, 0,self.startTime)
        deltas = np.diff( times )

        if not self.conditional:
            path += np.random.randn( N, len(times)-1 )*self.sigma*np.sqrt(deltas) + self.mu*deltas
            return path.cumsum(axis=1)
        
                
        else:
            """
            Alternatively, this can be accomplished by sampling directly from a multivariate normal given a linear
            projection. Ie
            
            N | N dot 1 = endPosition ~ Nor( 0, Sigma ), where Sigma is a diagonal matrix with elements proportional to
            the delta. This only problem with this is Sigma is too large for very large len(times).
            
            """
            T = self.endTime - self.startTime
            x = self.startTime
            for i, delta in enumerate( deltas ):
                x = x*(1-delta/T)+(self.endPosition - self.startPosition)*delta/T + self.sigma*np.sqrt(delta/T*(T-delta))*self.Nor.rvs(N)
                T = T - delta
                path[:,i] = x
            if abs(T -0)<1e-10:
                path[:,-1] = (self.endPosition - self.startPosition)
            path = path + self.startPosition
        return path
    
    def generate_max(self,t):
        pass
    
    def generate_min(self,t):
        pass
    
    def drift(t,x):
        return self.mu
    def diffusion(t,x):
        return self.sigma
    def diffusion_prime(t,x):
        return 0

        
        
        
#have a Vasicek model that has an alternative parameterizatiom but essentially just maps to OU_process

class OU_process(Diffusion_process):
    
    """
    The Orstein-Uhlenbeck process is a mean reverting model that is often used in finance to model interest rates.
    It is also known as a Vasicek model. It is defined by the SDE:
    
    dOU_t = theta*(mu-OU_t)*dt + sigma*dB_t$
    
    
    The model flucuates around its long term mean mu. mu is also a good starting Position of the process.
    
    There exists a solution to the SDE, see wikipedia.
    
    parameters: 
        theta: float > 0
        mu: float
        sigma: float > 0
        
          
    """
    def __init__(self, theta, mu, sigma, startTime = 0, startPosition = 0, endPosition = None, endTime = None ):
        assert sigma > 0 and theta > 0, "theta > 0 and sigma > 0."
        super(OU_process, self).__init__(startTime, startPosition, endTime, endPosition)
        self.theta = theta
        self.mu = mu
        self.sigma = sigma
        self.Normal = stats.norm()


    def _mean(self,t):
        if self.conditional:
            return super(OU_process,self)._mean(t) #TODO
        else:
            return self.startPosition*np.exp(-self.theta*(t-self.startTime))+self.mu*(1-np.exp(-self.theta*(t-self.startTime)))
                                                                
    def _var(self,t):
        if self.conditional:
            return super(OU_process,self)._get_variance_at(t)
        else:
            return self.sigma**2*(1-np.exp(-2*self.theta*t))/(2*self.theta)
    
    def _transition_pdf(self,x,t,y):
            mu = x*np.exp(-self.theta*t)+self.mu*(1-np.exp(-self.theta*t))
            sigmaSq = self.sigma**2*(1-np.exp(-self.theta*2*t))/(2*self.theta)
            return np.exp(-(y-mu)**2/(2*sigmaSq))/np.sqrt(2*pi*sigmaSq)            
    
    
    def _sample_position(self,t):
        if not self.conditional:
            return self.get_mean_at(t)+np.sqrt(self.get_variance_at(t))*self.Normal.rvs()
        else:
            #this needs to be completed
            return super(OU_process,self)._generate_position_at(t)
    
    def sample_path(self,times, N= 1, return_normals = False):
        "the parameter Normals = 0 is used for the Integrated OU Process"
        if not self.conditional:
            path=np.zeros( (N,len(times)) )
            times = np.insert( times, 0,self.startTime)
            path[ :, 0] = self.startPosition
            
            deltas = np.diff(times )
            normals = np.random.randn( N, len(times)-1 )
            x = self.startPosition*np.ones( N) 
            sigma = np.sqrt(self.sigma**2*(1-np.exp(-2*self.theta*deltas))/(2*self.theta))
            for i, delta in enumerate(deltas):
                mu = self.mu + np.exp(-self.theta*delta)*(x-self.mu)
                path[:, i] = mu + sigma[i]*normals[:,i] 
                x = path[:,i]
                """
                It would be really cool if there was a numpy func like np.cumf( func, array )
                that applies func(next_x, prev_x) to each element. For example, lambda x,y: y + x
                is the cumsum, and lambda x,y: x*y is the cumprod function.
                """
            if return_normals:
                return (path, normals )
            else:
                return path
        
        else:
            #TODO
            path = bridge_creation(self,times)
            return path
    
    def drift(self, x, t):
        return self.theta*(self.mu-x)
    def diffusion( self, x,t ):
        return self.sigma
        
class Integrated_OU_process(Diffusion_process):
    """
    The time-integrated Orstein-Uhlenbeck process
    $IOU_t = IOU_0 + \int_0^t OU_s ds$
    where $dOU_t = \theta*(\mu-OU_t)*dt + \sigma*dB_t, 
    OU_0 = x0$
    
    parameters:
    {theta:scalar > 0, mu:scalar, sigma:scalar>0, x0:scalar}
    
    modified from http://www.fisica.uniud.it/~milotti/DidatticaTS/Segnali/Gillespie_1996.pdf
    """
    def __init__(self, parameters, time_space_constraints):
        super(Integrated_OU_process,self).__init__( time_space_constraints)
        self.OU = OU_process({"theta":parameters["theta"], "mu":parameters["mu"], "sigma":parameters["sigma"]}, {"startTime":time_space_constraints["startTime"], "startPosition":parameters["x0"]})
        for p in parameters:
            setattr(self, p, float(parameters[p]))
        self.Normal = stats.norm()
        
    def _get_mean_at(self,t):
        delta = t - self.startTime
        if self.conditional:
            pass
        else:
            return self.startPosition + (self.x0-self.mu)/self.theta + self.mu*delta\
                            -(self.x0-self.mu)*np.exp(-self.theta*delta)/self.theta
    
    def _get_variance_at(self,t):
        delta = t - self.startTime
        if self.conditional:
            pass
        else:
            return self.sigma**2*(2*self.theta*delta-3+4*np.exp(-self.theta*delta)
                                  -2*np.exp(-2*self.theta*delta))/(2*self.sigma**3)


    def _generate_position_at(self,t):
        if self.conditional:
            pass 
        else:
            return self.get_mean_at(t)+np.sqrt(self.get_variance_at(t))*self.Normal.rvs()
    
    def _transition_pdf(self,x,t,y):
        mu = x + (self.x0 - self.mu)/self.theta + self.mu*t - (self.x0-self.mu)*np.exp(-self.theta*t)/self.theta
        sigmaSq = self.sigma**2*(2*self.theta*t-3+4*np.exp(-self.theta*t)-2*np.exp(-2*self.theta*t))/(2*self.sigma**3)
        return np.exp(-(y-mu)**2/(2*sigmaSq))/np.sqrt(2*pi*sigmaSq) 
    
        
    def generate_sample_path(self,times, returnUO = 0):
        "set returnUO to 1 to return the underlying UO path as well as the integrated UO path."
        if not self.conditional: 
            xPath, listOfNormals = self.OU.generate_sample_path(times, 1)
            path = []
            t = self.startTime
            y = self.startPosition
            for i, position in enumerate(xPath):
                delta = position[0]-t
                x = position[1]
                if delta != 0:
                    #there is an error here, I can smell it.
                    sigmaX = self.sigma**2*(1-np.exp(-2*self.theta*delta))/(2*self.theta)
                    sigmaY = self.sigma**2*(2*self.theta*delta-3+4*np.exp(-self.theta*delta)
                                      -np.exp(-2*self.theta*delta))/(2*self.sigma**3)
                    muY = y + (x-self.mu)/self.theta + self.mu*delta-(x-self.mu)*np.exp(-self.theta*delta)/self.theta
                    covXY = self.sigma**2*(1+np.exp(-2*self.theta*delta)-2*np.exp(-self.theta*delta))/(2*self.theta**2)
                    y = muY + np.sqrt(sigmaY - covXY**2/sigmaX)*self.Normal.rvs()+ covXY/np.sqrt(sigmaX)*listOfNormals[i]
                t = position[0]
                path.append((t,y))
            if returnUO==0:
                return path
            else:
                return path, xPath
        else:
           path = bridge_creation(self,times)
           if returnUO==0:
               return path
           else:
               return path, xPath 
        
    def _process2latex(self):

        return """$IOU_t = IOU_0 + \int_0^t OU_s ds
            \text{where} dOU_t = %.3f(%.3f-OU_t)dt + %.3fdB_t, 
            OU_0 = x0 
            """ %(self.theta, self.mu, self.sigma)
           
class SqBessel_process(Diffusion_process):
    """
    The (lambda0 dimensional) squared Bessel process is defined by the SDE:
    dX_t = lambda_0*dt + nu*sqrt(X_t)dB_t
    
    parameters:
    
        lambda_0: float,
        nu: float > 0
    
     Attn: startPosition and endPosition>0
     
     Based on R.N. Makarov and D. Glew's research on simulating squared bessel process. See "Exact
     Simulation of Bessel Diffusions", 2011.
     
    """
    def __init__(self, lambda_0, nu, startTime = 0, startPosition = 1, endTime = None, endPosition = None ):
            super(SqBessel_process, self).__init__(startTime, startPosition, endTime, endPosition)
            try:
                self.endPosition = 4.0/parameters["nu"]**2*self.endPosition
                self.x_T = self.endPosition
            except:
                pass
            self.x_0 = 4.0/nu**2*self.startPosition
            self.nu = nu
            self.lambda0 = lambda_0
            self.mu = 2*float(self.lambda0)/(self.nu*self.nu)-1
            self.Poi = stats.poisson
            self.Gamma = stats.gamma
            self.Nor = stats.norm
            self.InGamma = IncompleteGamma
            

    
    def generate_sample_path(self,times,absb=0):
        """
        absb is a boolean, true if absorbtion at 0, false else. See class' __doc__ for when 
        absorbtion is valid.
        """
        if absb:
            return self._generate_sample_path_with_absorption(times)
        else:
            return self._generate_sample_path_no_absorption(times)
       
    
    def _transition_pdf(self,x,t,y):
        try:
            return (y/x)**(0.5*self.mu)*np.exp(-0.5*(x+y)/self.nu**2/t)/(0.5*self.nu**2*t)*iv(abs(self.mu),4*np.sqrt(x*y)/(self.nu**2*t))
        except AttributeError:
            print("Attn: nu must be known and defined to calculate the transition pdf.")
            
    def _generate_sample_path_no_absorption(self, times):
        "mu must be greater than -1. The parameter times is a list of times to sample at."
        if self.mu<=-1:
            print("Attn: mu must be greater than -1. It is currently %f."%self.mu)
            return
        else:
            if not self.conditional:
                x=self.startPosition
                t=self.startTime
                path=[]
                for time in times:
                    delta=float(time-t)
                    try:
                        y=self.Poi.rvs(0.5*x/delta)
                        x=self.Gamma.rvs(y+self.mu+1)*2*delta
                    except:
                        pass
                    path.append((time,x))
                    t=time
            else:
                path = bridge_creation(self, times, 0)
                return path
            return [(p[0],self.rescalePath(p[1])) for p in path]
         
 
        
    def _generate_sample_path_with_absorption(self,times):
        "mu must be less than 0."
        if self.mu>=0:
            print("Attn: mu must be less than 0. It is currently %f."%self.mu)
        else:
            if not self.conditional:
                path=[]
                X=self.startPosition
                t=self.startTime
                tauEst=times[-1]+1
                for time in times:
                    delta = float(time - t)
                    if tauEst>times[-1]:
                        p_a = gammaincc(abs(self.mu),0.5*X/(delta))
                        if np.random.rand() < p_a:
                            tauEst = time
                    if time<tauEst:
                        Y = self.InGamma.rvs(abs(self.mu),0.5*X/(delta))
                        X = self.Gamma.rvs(Y+1)*2*delta
                    else:
                        X=0
                    t=time
                    path.append((t,X))
            else:
                path = bridge_creation(self, times, 1)
            return  [(p[0],self.rescalePath(p[1])) for p in path]
                    

    def _generate_position_at(self,t):
        p = self.generate_sample_path([t])
        return p[0][1]
    
       
    def generate_sample_FHT_bridge(self,times):
        "mu must be less than 0. This process has absorption at L=0. It simulates the absorption at 0 at some random time, tao, and creates a bridge process."
        if self.mu>0:
            print("mu must be less than 0. It is currently %f."%self.mu)
        else:
            X=self.startPosition
            t=self.t_0
            path=[]
            FHT=self.startPosition/(2*self.Gamma.rvs(abs(self.mu)))
            for time in times:
                if time<FHT:
                    d=(FHT-t)*(time-t)
                    Y=self.Poi.rvs(X*(FHT-time)/(2*d))
                    X=self.Gamma.rvs(Y-self.mu+1)*d/(FHT-t)
                else:
                    X=0
                t=time
                path.append((t,X))
            return  [(p[0],self.rescalePath(p[1])) for p in path]
    
    
    def rescalePath(self,x):
        #All of the simulation algorithms assume nu=2, so we must
        # rescale the process to output a path that is has the user specified
        # nu. Note that this rescaling does not change mu.
        return self.nu**2/4.0*x
    
    
class CIR_process(Diffusion_process):
    """
    The CIR process is defined by
    dCIR_t = (lambda_0 - lambda_1*CIR_t)dt + nu*np.sqrt(CIR_t)*dB_t
    
    Due to the nature of this interface, absorption at 0 is impossible. See the PyProcess library for 
    the ability to aborb at 0.
    This is a mean reverting process if both lambdas are positive: the process flucuates around lambda_0/lambda_1
    
    parameters:
    {lambda_0:scalar, lambda_1:scalar, nu:scalar>0}
    
    """
    def __init__(self, parameters, space_time_constraints):
        super(CIR_process,self).__init__(space_time_constraints)
        for p in parameters:
            setattr(self, p, float(parameters[p]))
        self.Normal = stats.norm()
        #transform the space time positions 
        _space_time_constraints = {}
        _space_time_constraints['startTime'] = self._time_transformation(space_time_constraints['startTime'])
        _space_time_constraints['startPosition'] = self._inverse_space_transformation(space_time_constraints['startTime'], space_time_constraints['startPosition'])
        try: 
            _space_time_constraints['endPosition'] = self._inverse_space_transformation(space_time_constraints['endTime'], space_time_constraints['endPosition'])
            _space_time_constraints['endTime'] = self._time_transformation(space_time_constraints['endTime'])
        except:
            pass
        self.SqB = SqBessel_process({"lambda0":parameters["lambda_0"], "nu":parameters["nu"]}, _space_time_constraints) #need to change start position for non-zero startTime
        self.mu=self.SqB.mu
    
    def _process2latex(self):
        return """
        $dCIR_t = (%.3f - %.3fCIR_t)dt + %.3f\np.sqrt(CIR_t)dB_t$
        """%(self.lambda_0, self.lambda_1, self.nu)
        
    def _transition_pdf(self,x,t,y):
        return np.exp(self.lambda_1*t)*SqB._transition_pdf(x,self._time_transformation(t), np.exp(self.lambda_1*t)*y)
    
    def generate_sample_path(self, times, abs=0):
        "abs is a boolean: true if desire nonzero probability of absorption at 0, false else."
        #first, transform times:
        transformedTimes = [self._time_transformation(t) for t in times]
        path = self.SqB.generate_sample_path(transformedTimes,abs)
        tpath = [self._space_transformation(times[i],p[1]) for i,p in enumerate(path) ]
        path=[]
        for i in xrange(len(tpath)):
            path.append((times[i],tpath[i]))
        return path
    
    def _generate_position_at(self,t):
        t_prime = self._time_transformation(t)
        x = self.SqB.generate_position_at(t_prime)
        return self._space_transformation(t,x)
    
    def _time_transformation(self,t):
        if self.lambda_1==0:
            return t
        else:
            return (np.exp(self.lambda_1*t)-1)/self.lambda_1
        
    def _space_transformation(self,t,x):
        return np.exp(-self.lambda_1*t)*x
    
    def _inverse_space_transformation(self,t,x):
        return np.exp(self.lambda_1*t)*x
        
    def _inverse_time_transformation(self, t):
        if self.lambda_1==0:
            return t
        else:
            return np.log(self.lambda_1*t+1)/self.lambda_1
            
    def _get_mean_at(self,t):
        pass
    
    def _get_variance_at(self,t):
        pass

class CEV_process(Diffusion_process):
    """
    defined by:
    $$dCEV = rCEVdt + \deltaCEV^{\beta+1}dW_t$$
    
    parameters:
    {r: scalar, delta:scalar>0, beta:scalar<0} #typically beta<=-1/2
    
    """
    
    def __init__(self,parameters, time_space_constraints):
        super(CEV_process,self).__init__(time_space_constraints)
        for p in parameters:
            setattr(self, p, parameters[p])
        time_space_constraints["startPosition"]=self.CEV_to_SqB(self.startPosition)
        self.SqB = SqBessel_process({"lambda0":(2-1/self.beta), "nu":2}, time_space_constraints)
        
    def _process2latex(self):
        return 
        """
        $dCEV_t = %.3fCEVdt + %.3fCEV^{\%.3f + 1}dB_t$
        """%(self.r, self.delta, self.beta)
        
    def _time_transform(self,t):
        if self.r*self.beta==0:
            return t
        else:
            return (np.exp(self.r*self.beta*2*t)-1)/(self.r*self.beta*2)
        
    def CEV_to_SqB(self,x):
        return x**(-2*self.beta)/(self.delta*self.beta)**2
    
    def _scalar_space_transform(self,t,x):
        return np.exp(self.r*t)*x
    
    def SqB_to_CEV(self,x):
        ans = (self.delta**2*self.beta**2*x)**(-1/(2.0*self.beta))
        return ans
    
    
    def generate_sample_path(self,times, abs=0):
        if self.r==0:
            SqBPath = self.SqB.generate_sample_path(times, abs)
            return [(x[0],self.SqB_to_CEV(x[1])) for x in SqBPath]
        else:
            transformedTimes = [self._time_transform(t) for t in times]
            SqBPath = self.SqB.generate_sample_path(transformedTimes, abs)
            tempPath = [self.SqB_to_CEV(x[1]) for x in SqBPath]
            return [(times[i], self._scalar_space_transform(times[i], p) ) for i,p in enumerate(tempPath)]
        
    def _generate_position_at(self,t):
        if self.r==0:
            SqBpos = self.SqB.generate_position_at(t)
            return self.SqB_to_CEV(SqBpos)
        else:
            transformedTime = self._time_transform(t)
            SqBpos = self.SqB.generate_position_at(transformedTime)
            return self._scalar_space_transform(t,self.SqB_to_CEV(SqBpos))
         
class Periodic_drift_process(Diffusion_process):
    """
      dX_t = psi*sin(X_t + theta)dt + dBt 
      
      parameters:
      {psi:scalar>0, theta:scalar>0}
      
      
      This cannot be conditioned on start or end conditions.
      Extensions to come. 
    """
    def __init__(self, parameters, space_time_constraints):
        """Note that space-time constraints cannot be given"""
        space_time_constraints = {"startTime":0, "startPosition":0}
        super(Periodic_drift_process,self).__init__(space_time_constraints)
        self.psi = parameters["psi"]
        self.theta = parameters["theta"]
        self._findBounds()        
        self.BB = Wiener_process({"mu":0, "sigma":1}, space_time_constraints)
        self.Poi = Marked_poisson_process({"rate":1, "U":self.max, "L":self.min, "startTime":0}) #need to create marked poisson process class
        self.Nor = stats.norm
        self.Uni = stats.uniform()

    def _process2latex(self):
        return """
        $dX_t = %.3f*sin(X_t + %.3f)dt + dB_t$
        """%(self.psi, self.theta)
    
    

    def __generate_sample_path(self,T,x):
        #generates a path of length 2. This is for efficiency issues
        self.BB.startPosition=x
        while (True):
            
            
            #produce a marked poisson process
            markedProcess = self.Poi.generate_marked_process(T)

            #generate end point using an AR scheme
            while (True):
                N = x+ np.sqrt(T)*self.Nor.rvs()
                U = self.Uni.rvs()
                if U<=np.exp(-self.psi*cos(N-self.theta)+self.psi):
                    break
            self.BB.endPosition = N
        
            #generate brownian bridge
            try:
                skeleton = self.BB.generate_sample_path([p[0] for p in markedProcess])
            except:
                skeleton = []
            transformSkeleton = [self.phi(p[1]) for p in skeleton]
            
            #calculate indicators
            I=1
            for i in xrange(len(transformSkeleton)):
                if transformSkeleton[i]>markedProcess[i][1]:
                    I=0
                    break
            
            #check indicators
            if I==1:
                return N, skeleton + [(T,N)]
                
    def generate_sample_path(self,times):
        "currently will only return a random path before time T. Can be connected by brownian bridges"
        #this algorithm uses the EA1 algorithm by Beskos and 
        # Roberts on exact simulation of diffusions. It's an AR algorithm.
        # For some parameters, the probability of acceptance can be very low.
        time = 0
        endPoint = 0
        skeleton=[]
        T = times[-1]
        while time<T:
            if time+2<T:
                delta=2
            else:
                delta = T-time
            endPoint, tempSkeleton = self.__generate_sample_path(delta,endPoint)
            tk = [(time + x[0],x[1]) for x in tempSkeleton] 
            skeleton+=tk
            
            time+=2
        
     
        
        return self._construct_path_from_skeleton([(0,self.startPosition)]+skeleton, times)

    def _construct_path_from_skeleton(self, skeleton, times):
        i=0
        path=[]
        self.BB.startPosition = self.startPosition
        self.BB.startTime = self.startTime
        for (tao,x) in skeleton[1:]:
            self.BB.endPosition = x
            self.BB.endTime = tao
            temptimes = []
            while ( i<len(times) ) and (times[i] <= tao):
                temptimes.append(times[i])
                i+=1
            for p in self.BB.generate_sample_path(temptimes):
                path.append(p)
            self.BB.startPosition = x
            self.BB.startTime = tao
        return path
        
    
    def _findBounds(self):
        if self.psi<=.5:
            self.max = self.psi
        else:
            self.max = 0.125+0.5*(self.psi*self.psi+self.psi)
        self.min = -self.psi*.5
    
    def phi(self,x):
        return (0.5*self.psi*self.psi*sin(x-self.theta)**2+0.5*self.psi*cos(x-self.theta)-self.min)/self.max

class GBM_process(Diffusion_process):
    """
    Geometric Brownian Motion process is defined by the SDE 
    dGBM_t = mu GBM_tdt + sigma GBM_t dW_t
    
    and has general solution 
    GBM_t = GBM_0 exp( (mu - 0.5sigma^2)t + sigma W_t) 
    You can use this to reparamertize the sde above.
    
    parameters:
        mu: float
        sigma: float
        
    """
    
    def __init__(self, mu, sigma, startPosition = 1, startTime = 0, endPosition = None, endTime = None):
        super(GBM_process,self).__init__(startTime, startPosition, endTime, endPosition)
        self.mu = mu
        self.sigma = sigma
        self.Nor = stats.norm(0,1)
        try:
            wienerConstraints = {"startTime":self.startTime, "endTime":self.endTime, "startPosition":0, "endPosition":np.log(self.endPosition/self.startPosition)}
        except:
            wienerConstraints = {"startTime":self.startTime, "startPosition":0}
        self.wiener = Wiener_process( mu = (self.mu-0.5*self.sigma**2), sigma = self.sigma, **wienerConstraints)
        

    
    def _mean(self,t):
        if not self.conditional:
            return self.startPosition*np.exp(self.mu*t)
        else:
           delta = self.endPosition - self.startPosition
           return self.startPosition*np.exp(-0.5*(self.sigma*t)**2/delta+self.Weiner.endPosition*t/delta)        
            
    def _transition_pdf(self,x,t,y):
        delta = t - self.startTimee
        d = (self.mu - 0.5*self.sigma**2)
        return  x/(y*np.sqrt(2*Pi*self.sigma**2**delta))*np.exp(-(np.log(y/x)-d)**2/(2*self.sigma**2*delta))
    
    def _var(self,t):
        if not self.conditional:
            return self.startPosition**2*np.exp(2*self.mu*t)*(np.exp(self.sigma**2*t)-1)
        else:
            X=np.log(self.endPosition/self.startPosition)-(self.mu-0.5*(self.sigma**2))*T
            X=X/self.sigma
            delta = self.endTime-self.startTime
            return self.startPosition**2*np.exp(2*self.mu*t-self.sigma**2*t + self.sigma*t/delta*(self.sigma*(delta-T)+2*X))*(np.exp(self.sigma**2*(delta-t)*t/delta)-1)
        
    def _sample_position(self,t, N = 1):
        return self.startPosition*np.exp( self.wiener.sample_position(t,N) )
        
    def _sample_path(self,times, N = 1):
        return self.startPosition*np.exp(self.wiener.sample_path(times, N) )
        
    def drift(self, x, t):
        return self.mu*x
        
    def diffusion(self,x,t):
        return self.sigma*x
        
        
class Custom_diffusion(Diffusion_process):
    """
    simulates the diffusion:
    dX_t = a(X_t,t)*dt + b(X_t,t)*dW_t
    
    where a and b are inputed functions
    
    parameters:
    {a:time and space function, b:nonnegative time and space function}
    Ex:
    def f1(x,t):
        return x^2-t
        
    def f2(x,t):
        return np.sqrt(x)
    param = {"a":f1, "b":f2}
    
    
    
    The superclass Diffusion_process contains most of the methods for this class. See it's documentation for all.
    """
    def __init__(self,parameters,space_time_constraints):
        super(Custom_diffusion, self).__init__(space_time_constraints)
        self.drift = parameters["a"]
        self.diffusion = parameters["b"]
    
    
    
    

class Jump_Diffusion_process(object):
    
    def __init__(self,dict):    
        try:
            self.startTime = dict["startTime"]
            self.startPosition = dict["startPosition"]
            self.conditional=False
            if dict.has_key("endTime"):
                self.endTime = dict["endTime"]
                self.endPosition = dict["endPosition"]
                self.conditional=True
        except KeyError:
            print("Missing constraint in initial\end value dictionary. Check spelling?" )
    def transition_pdf(self,t,y):
        "this method calls self._transition_pdf(x) in the subclass"
        self._check_time(t)
        if not self.conditional:
            return self._transition_pdf(self.startPosition, t-self.startTime, y)
        else:
            return self._transition_pdf(self.startPosition, t-self.startTime, y)*self._transition_pdf(y, self.endTime-t, self.endPosition)\
                    /self._transition_pdf(self.startPosition,self.endTime - self.startTime, self.endPosition)
    

    def expected_value(self,f,t,N):
        "uses a monte carlo approach to evaluate the expected value of the process f(X_t). N is the number of iterations. The parameter f\
        is a univariate python function."
        print("Attn: performing a Monte Carlo simulation...")
        self._check_time(t)
        if not self.conditional:
            sum=0
            for i in xrange(N):
                sum+=f(self.generate_position_at(t))
            return sum/N
        else:
            #This uses a change of measure technique.
            sum=0
            self.conditional=False
            for i in xrange(N):
                X = self.generate_position_at(t)
                sum+=self._transition_pdf(X,self.endTime-t,self.endPosition)*f(X)
            self.conditional=True
            return sum/(N*self._transition_pdf(self.startPosition, self.endTime-self.startTime, self.endPosition))
        
    def generate_position_at(self,t):
        "if _get_position_at() is not overwritten in a subclass, this function will use euler scheme"
        self._check_time(t)        
        if self.startTime<t:
            return self._generate_position_at(t)
        
    
    def get_mean_at(self,t):
        self._check_time(t)
        return self._get_mean_at(t)
    
    
    def get_variance_at(self,t):
        self._check_time(t)
        return self._get_variance_at(t)
    
    def generate_sample_path(self,times):
        '"times" is a list of times, with the first time greater than initial time."'
        self._check_time(times[0])
        return self._generate_sample_path(times)
    
    
    def _generate_sample_path(self,times):
        "I need some way to use a euler scheme AND evaluate the approximation at t in times"
        pass
    
    def _get_variance_at(self,t):
        var = SampleVarStat()
        simulations = 10000
        for i in xrange(simulations):
            var.push(self.generate_position_at(t))
        return var.get_variance_at()
    
    
    def _get_mean_at(self,t):
        "if _get_mean_at() is not overwritten, then we use MC methods; 100000 iterations"
        def id(x):
            return x
        return self.expected_value(id, t, 100000)
    
    def _generate_position_at(self,T,delta=0.001):
        return self.Euler_scheme(t,return_array=False)
        
    def _transition_pdf(self,x,t,y):
        print ("Attn: transition pdf not defined")

    def _check_time(self,t):
        if t<self.startTime:
            print ("Attn: inputed time not valid (check if beginning time is less than startTime).")
    
        

class Gamma_process(Jump_Diffusion_process):
    """
    Defined by G(t+h) - G(t) is distributed as a Gamma random variable. 
    parameters:
    {"mean":scalar>0, "variance":scalar>0}
    or
    {"rate":scalar>0, "size":scalar>0}
    
    
    Under the first parameterization: E[G(t)] = mean*t, Var(G(t)) = variance*t
    Under the second parameterization: G(t+1) - G(t) ~ Gamma(rate, size)
    where Gamma has pdf [size^(-rate)/gamma(rate)]*x^(rate-1)*np.exp(-x/size)
    
    see http://eprints.soton.ac.uk/55793/1/wsc03vg.pdf for details on this process.
    
    """
    
    def __init__(self, parameters, time_space_constraints):
        super(Gamma_process,self).__init__(time_space_constraints)
        try:
            self.mean = float(parameters["rate"]/parameters["size"])
            self.variance = float(self.mean/parameters["size"])
        except KeyError:
            self.mean = float(parameters["mean"])
            self.variance = float(parameters["variance"])
        self.gamma = stats.gamma
        self.beta = stats.beta
        
    def _generate_position_at(self,t):
        if not self.conditional:
            return self.startPosition + self.gamma.rvs(self.mean**2*(t-self.startTime)/self.variance)*self.mean/self.variance
        else:
            return self.startPosition + (self.endPosition-self.startPosition)*self.beta.rvs((t-self.startTime)/self.variance, (self.endTime-t)/self.variance)
            
    def _get_mean_at(self,t):
        "notice the conditional is independent of self.mean."
        if self.conditional:
            return self.startPosition + (self.endPosition-self.startPosition)*(t- self.startTime)/(self.endTime - self.startTime)
        else:
            return self.startPosition + self.mean*(t-self.startTime)
    
    def _get_variance_at(self,t):
        if self.conditional:
            alpha = (t-self.startTime)/self.variance
            beta = (self.endTime-t)/self.variance
            return (self.endPosition-self.startPosition)**2(alpha*beta)/(alpha+beta)**2/(alpha+beta+1)
        else:
            return self.variance*(t-self.startTime)
    
        
    def _generate_sample_path(self,times):
        if not self.conditional:
            t = self.startTime
            x = self.startPosition
            path=[]
            for time in times:
                delta = time - t
                try:
                    g = self.gamma.rvs(self.mean**2*delta/self.variance)*self.variance/self.mean
                    x = x + g
                except ValueError:
                    pass
                t = time
                path.append((t,x))
            return path
        else:
            x = self.startPosition
            t = self.startTime
            path=[]
            for time in times:
                delta1 = time -t
                delta2 = self.endTime - time
                if (delta1!=0 and delta2!=0):
                    b = (self.endPosition-x)*self.beta.rvs(delta1/self.variance,delta2/self.variance)
                elif delta1 == 0:
                    b = 0
                else:
                    b = self.endTime - x    
                x = x + b
                t = time
                path.append((t,x))
            return path
    
    def _transition_pdf(self,x,t,y):
        return self.variance/self.mean*self.gamma.pdf((y-x)*self.mean/self.variance*t,self.mean**2/self.variance)
    
    
class Gamma_variance_process(Jump_Diffusion_process):
    """
    The Gamma variance process is a brownian motion subordinator:
    VG_t = mu*G_t(t,a,b) + sigma*B_{G_t(t,a,b)}
    i.e. VG process is a time transformed brownian motion plus a scalar drift term.
    It can also be represented by the difference of two gamma processes. 
    
    Note: currently, if conditional, endPosition=0. Further extensions will be to make this nonzero.
    
    parameters:
    {mu:scalar, sigma:scalar>0, variance:scalar>0} 
    or
    {mu:scalar, sigma:scalar>0, rate:scalar>0}
    *note the mean of the gamma process is 1, hence the reduction of parameters needed.
    
    The parameterization depends on the way you parameterize the underlying gamma process.   
    See http://www.math.nyu.edu/research/carrp/papers/pdf/VGEFRpub.pdf"
    """
    
    def __init__(self,parameters, space_time_constraints):
        super(Gamma_variance_process,self).__init__(space_time_constraints)
        try:
            self.variance=v = float(parameters["variance"])
        except KeyError:
            self.variance=v = 1/float(parameters["rate"])
        self.sigma=s = float(parameters["sigma"])
        self.mu= m = float(parameters["mu"])
        self.mu1 = 0.5*np.sqrt(m**2 + 2*s**2/v)+m/2
        self.mu2 = 0.5*np.sqrt(m**2 + 2*s**2/v)-m/2
        self.var1 = self.mu1**2*v
        self.var2 = self.mu2**2*v
        self.GamPos = Gamma_process({"mean":self.mu1,"variance":self.var1}, space_time_constraints)
        self.GamNeg = Gamma_process({"mean":self.mu2, "variance":self.var2}, space_time_constraints)
        
        
    def _get_mean_at(self,t):
        return self.GamPos.get_mean_at(t)-self.GamNeg.get_mean_at(t)
    
    def _get_variance_at(self,t):
        return self.GamPos.get_variance_at(t)+self.GamNeg.get_variance_at(t)
        
    def _generate_position_at(self,t):
        x = self.GamPos.generate_position_at(t)
        y = self.GamNeg.generate_position_at(t)
        return x-y

    def _generate_sample_path(self,times):
        path=[]
        pathPos = self.GamPos.generate_sample_path(times)
        pathNeg = self.GamNeg.generate_sample_path(times)
        for i,time in enumerate(times):
            path.append((time, pathPos[i][1] - pathNeg[i][1]))
        return path
    
    def _transition_pdf(self,x,t,y):
        alpha1=self.mu1**2*t/self.var1
        beta1 = self.var1/self.mu2
        alpha2 = self.mu2**2*t/self.var2
        beta2 = self.var2/self.mu2
        return np.exp(-(y-x))
    
    
class Geometric_gamma_process(Jump_Diffusion_process):
        """
        the geometric gamma process has the representation GG_0*np.exp(G_t) where G_t is a gamma process.
        
        parameters
        {mu: scalar, sigma: scalar>0 } 
        The parameters refer to the parameters in the gamma process (mu = mean, sigma = variance) (see gamma_process documentation for more details).
    
        """
        
        def __init__(self, parameters, space_time_constraints):
            super(Geometric_gamma_process, self).__init__(space_time_constraints)
            self.mu = float(parameters["mu"])
            self.sigma = float(parameters["sigma"])
            try:
                self.gammaProcess = Gamma_process({"mean":self.mu, "variance":self.sigma}, {"startPosition":0, "startTime":0, "endTime":self.endTime, "endPosition":np.log(self.endPosition/self.startPosition)})
            except:
                self.gammaProcess = Gamma_process({"mean":self.mu, "variance":self.sigma}, {"startPosition":0, "startTime":0})

        
        
        def _transition_pdf(self,z,t,x):
            "as this is a strictly increasing process, the condition x<z must hold"
            alpha = self.mu**2*t/self.sigma
            beta = self.sigma/self.mu
            return (np.log((z-x)/self.startPosition))**(alpha-1)*((z-x)/self.startPosition)**(-1/beta)/((z-x)*beta**alpha*gamma(alpha))
        
        def _get_mean_at(self,t):
            if self.conditional:
                pass
            else:
                if self.sigma<self.mu: #this condition gauruntees the expectation exists
                    return self.startPosition*(1-self.sigma/self.mu)**(-self.mu**2*(t-self.startTime)/self.sigma)
                else:
                    print ("Attn: does not exist.")
                    return "DNE"
        
        def _generate_sample_path(self,times):
            return [(p[0],self.startPosition*np.exp(p[1])) for p in self.gammaProcess.generate_sample_path(times)]
        
        def _generate_position_at(self,t):
            return np.exp(self.gammaProcess.generate_position_at(t))
        

class Inverse_Gaussian_process(Jump_Diffusion_process):
    """
    Based on 
    http://finance.math.ucalgary.ca/papers/CliffTalk26March09.pdf
    params:
    {a: >0  , b: >0 }
    
    Note: currently cannot be conditioned on endTime nor endPosition.
    """
    
    def __init__(self, parameters, timespace_constraints):
        super(Inverse_Gaussian_process,self).__init__(timespace_constraints)
        self.a = parameters['a']
        self.b = parameters['b']
        self.IG = InverseGaussian()
        
        
    def _generate_position_at(self, t):
        return self.startPosition + self.IG.rvs(self.a*(t - self.startTime),self.b) 
    
    def _generate_sample_path(self, times):
        x = self.startPosition
        t = self.startTime
        path = []
        for time in times:
                delta = time - t
                x += self.IG.rvs(self.a*delta, self.b)
                t = time
                path.append((t,x))
        return path
    
    def _get_mean_at(self,t):
        return self.startPosition + self.a*t/self.b
    
    def _get_variance_at(self,t):
        return (self.a*t/self.b)**3/(self.a*t)
    
    def _transition_pdf(self,z,t,x):
        y = x-z
        return self.a*(t-self.startTime)/(2*pi)*(y)**(-1.5)*np.np.exp(-0.5*( (self.a*(t-self.startTime))**2/y + self.b**2*y ) + self.a*(t-self.startTime)*self.b )


class Normal_Inverse_Gaussian_process(Jump_Diffusion_process):
    """This is a Brownian motion subordinated by a inverse gaussian process.
    From http://finance.math.ucalgary.ca/papers/CliffTalk26March09.pdf
    
    paramters 
    {beta: scalar, alpha: |beta|<alpha, delta:>0}
    """
    
    def __init__(self, parameters, timespace_constraints):
        super(Normal_Inverse_Gaussian_process,self).__init__(timespace_constraints)
        self.alpha = parameters['alpha']
        self.beta = parameters['beta']
        self.delta = parameters['delta']
        
        b = self.delta*np.sqrt(self.alpha**2 - self.beta**2)
        self.IGprocess = Inverse_Gaussian_process( {"a":1, "b":b}, {"startTime":0, "startPosition":self.startTime})
        self.BMprocess = Wiener_process({'mu':self.beta*self.delta**2, 'sigma':self.delta}, 
                                        {"startTime":self.startTime, 'startPosition':self.startPosition })
    
    def _generate_position_at(self, t):
        return self.BMprocess.generate_position_at(self.IGprocess.generate_position_at(t))
    
    def _generate_sample_path(self, times):
        p = self.IGprocess.generate_sample_path(times)
        return self.BMprocess.generate_sample_path([x[1] for x in p])
    
    def _get_mean_at(self,t):
        return self.startPosition + self.delta*(t-self.startTime)*self.beta/(np.sqrt(self.alpha**2 - self.beta**2))
    
    def _get_variance_at(self,t):
        return self.delta*(t-self.startTime)*self.alpha**2/(np.sqrt(self.alpha**2 - self.beta**2))**3
    
    
    def _transition_pdf(self,z,t,x):
        y = np.sqrt(self.delta**2 - + (x-z)**2)
        gamma = np.sqrt(self.alpha**2 - self.beta**2)
        return self.alpha*self.delta*(t - self.startTime)*kn(1,self.alpha*y )*np.np.exp(self.delta*gamma+self.beta*(z-x) )/(pi*y)
    
    
    
    

class Custom_process(object):
    """
    This class is a user defined sum of processes. The parameters are classes from this module.    
    Ex: 
    WP = Wiener_process{parametersWP, spaceTimeConstraintsWP}
    PP = Poisson_process{parametersPP, space_time_constraintsPP}
    Custom = Custom_process( WP, PP )
    
    Custom.get_mean_at(10)
    
    
    """
    def __init__(self, *args):
        self.processes = []
        for arg in args:
            self.processes.append(arg)
        
    def get_mean_at(self,t):
        sum=0
        for p in self.processes:
            sum+=p.get_mean_at(t)
        return sum
    
    def get_variance_at(self,t):
        sum=0
        for p in self.processes:
            sum+=p.get_variance_at(t)
        return sum

    def generate_position_at(self,t):
        sum=0
        for p in self.processes:
            sum+=p.generate_position_at(t)
        return sum
    
    def generate_sample_path(self,times, *args):
        "returns an array of a path, not an immutable list!"
        path = [[t,0] for t in times]
        for i in xrange(len(self.processes)):
            try:
              tempPath = self.processes[i].generate_sample_path(times, args[i] )
            except:
              tempPath = self.processes[i].generate_sample_path(times)
            for k in xrange(len(tempPath)):
                path[k][1]+=tempPath[k][1]
        return path
    
    
    
    
            
            
#------------------------------------------------------------------------------ 
# auxilary classes and functions

class Constant(object):
    def __init__(self,c):
        self.c=c
        
    def mean(self):
        return self.c
    def var(self):
        return 0
    def rvs(self, n):
        return self.c*np.ones(n)
        
    def cdf(self,x):
        if x<self.c:
            return 0
        else:
            return 1
    def sf(self,x):
        return 1-self.cdf(x)
    
    def pmf(self,x):
        if x==self.c:
            return 1
        else:
            return 0


class InverseGaussian(object):
        
    def rvs(self, a, b):
        y = stats.norm.rvs()
        x = a/b + y/(2*b**2) + np.sqrt(4*a*b+ y**2)/(2*b**2)
        u = stats.uniform.rvs()
        if u <= a/(a+x*b):
            return x
        return (a/b)**2/x

def bridge_creation(process,times, *args):
   # this algorithm relies on the fact that 1-dim diffusion are time reversible.
    print ("Attn: using an AR method...")
    process.conditional = False
    temp = process.startPosition
    while (True):
        sample_path=[]
        forward = process.generate_sample_path(times, *args)
        process.startPosition = process.endPosition
        backward = process.generate_sample_path(reverse_times(process, times), *args)
        process.startPosition = temp
        check = (forward[0][1]-backward[-1][1]>0)
        i=1
        N = len(times)
        sample_path.append(forward[0])
        while (i<N-1) and (check == (forward[i][1]-backward[-1-i][1]>0) ):
            sample_path.append(forward[i])
            i+=1
        
        if i != N-1: #an intersection was found
            k=0
            while(N-1-i-k>=0):
                sample_path.append((times[i+k],backward[N-1-i-k][1]))
                k+=1
            process.conditional = True

            return sample_path
                
        
    
def reverse_times(process, times):
    reverse_times=[]
    for time in reversed(times):
        reverse_times.append(process.endTime - time - process.startTime)
    return reverse_times

def transform_path(path,f):
    "accepts a path, ie [(t,x_t)], and transforms it into [(t,f(x)]."
    return [(p[0],f(p[1])) for p in path]

class SampleVarStat(object):
    def __init__(self):
        self.S=0
        self.oldM=0
        self.newM=0
        self.k=1
        
        
    
    def push(self,x):
        if self.k==0:
            self.S=0
            self.oldM=x
        else:
            self.newM = self.oldM+(x-self.oldM)/self.k
            self.S+=(x-self.oldM)*(x-self.newM)
            self.oldM=self.newM
            
        self.k+=1
    
    def get_variance(self):
        return self.S/(self.k-1)
    
class IncompleteGamma(object):
    "defined on negative integers. This is untested for accuracy"
    "Used in Bessel process simulation."
    def __init__(self,shape=None,scale=None):
        self.shape = float(shape)
        self.scale = float(scale)
        self.Uni = stats.uniform
        
    def pdf(self,n,shape,scale):
        inc = gammainc(shape,scale) #incomplete gamma function
        return np.np.exp(-scale)*sp.power(scale, n+shape)/(gamma(n+1+shape)*inc)
    
    def cdf(self,n, shape,scale):
        sum=0
        for i in xrange(n+1):
            sum+=self.pdf(i,shape,scale)
        return sum
    
    def rvs(self,shape,scale):
        "uses a inversion method: the chop-down-search starting \
        at the mode"
        pos =  mode = float(max(0,int(shape-scale)))
        U = self.Uni.rvs()
        sum = self.pdf(mode,shape,scale)
        ub = mode+1
        lb = mode-1
        Pub = sum*scale/(mode+1+shape)
        Plb = sum*(mode+shape)/scale
        while sum<U:
            if Plb>Pub and lb>=0:
                sum+=Plb
                pos = lb
                Plb = (lb+shape)/scale*Plb
                lb-=1
            else:
                sum+=Pub
                pos= ub
                Pub = scale/(ub+1+shape)*Pub
                ub+=1
        return float(pos)
    
    def rvsII(self,shape,scale):
        U = self.Uni.rvs()
        pos = 0
        sum = self.pdf(pos,shape,scale)
        while sum<U:
            pos+=1
            sum+=self.pdf(pos,shape,scale)
        return float(pos)