'''getting started with diffusions, continuous time stochastic processes Author: josef-pktd License: BSD References ---------- An Algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations Author(s): Desmond J. Higham Source: SIAM Review, Vol. 43, No. 3 (Sep., 2001), pp. 525-546 Published by: Society for Industrial and Applied Mathematics Stable URL: http://www.jstor.org/stable/3649798 http://www.sitmo.com/ especially the formula collection Notes ----- OU process: use same trick for ARMA with constant (non-zero mean) and drift some of the processes have easy multivariate extensions *Open Issues* include xzero in returned sample or not? currently not *TODOS* * Milstein from Higham paper, for which processes does it apply * Maximum Likelihood estimation * more statistical properties (useful for tests) * helper functions for display and MonteCarlo summaries (also for testing/checking) * more processes for the menagerie (e.g. from empirical papers) * characteristic functions * transformations, non-linear e.g. log * special estimators, e.g. Ait Sahalia, empirical characteristic functions * fft examples * check naming of methods, "simulate", "sample", "simexact", ... ? stochastic volatility models: estimation unclear finance applications ? option pricing, interest rate models ''' import numpy as np from scipy import stats, signal import matplotlib.pyplot as plt #np.random.seed(987656789) class Diffusion: '''Wiener Process, Brownian Motion with mu=0 and sigma=1 ''' def __init__(self): pass def simulateW(self, nobs=100, T=1, dt=None, nrepl=1): '''generate sample of Wiener Process ''' dt = T*1.0/nobs t = np.linspace(dt, 1, nobs) dW = np.sqrt(dt)*np.random.normal(size=(nrepl, nobs)) W = np.cumsum(dW,1) self.dW = dW return W, t def expectedsim(self, func, nobs=100, T=1, dt=None, nrepl=1): '''get expectation of a function of a Wiener Process by simulation initially test example from ''' W, t = self.simulateW(nobs=nobs, T=T, dt=dt, nrepl=nrepl) U = func(t, W) Umean = U.mean(0) return U, Umean, t class AffineDiffusion(Diffusion): r''' differential equation: :math:: dx_t = f(t,x)dt + \sigma(t,x)dW_t integral: :math:: x_T = x_0 + \int_{0}^{T}f(t,S)dt + \int_0^T \sigma(t,S)dW_t TODO: check definition, affine, what about jump diffusion? ''' def __init__(self): pass def sim(self, nobs=100, T=1, dt=None, nrepl=1): # this does not look correct if drift or sig depend on x # see arithmetic BM W, t = self.simulateW(nobs=nobs, T=T, dt=dt, nrepl=nrepl) dx = self._drift() + self._sig() * W x = np.cumsum(dx,1) xmean = x.mean(0) return x, xmean, t def simEM(self, xzero=None, nobs=100, T=1, dt=None, nrepl=1, Tratio=4): ''' from Higham 2001 TODO: reverse parameterization to start with final nobs and DT TODO: check if I can skip the loop using my way from exactprocess problem might be Winc (reshape into 3d and sum) TODO: (later) check memory efficiency for large simulations ''' #TODO: reverse parameterization to start with final nobs and DT nobs = nobs * Tratio # simple way to change parameter # maybe wrong parameterization, # drift too large, variance too small ? which dt/Dt # _drift, _sig independent of dt is wrong if xzero is None: xzero = self.xzero if dt is None: dt = T*1.0/nobs W, t = self.simulateW(nobs=nobs, T=T, dt=dt, nrepl=nrepl) dW = self.dW t = np.linspace(dt, 1, nobs) Dt = Tratio*dt L = nobs/Tratio # L EM steps of size Dt = R*dt Xem = np.zeros((nrepl,L)) # preallocate for efficiency Xtemp = xzero Xem[:,0] = xzero for j in np.arange(1,L): #Winc = np.sum(dW[:,Tratio*(j-1)+1:Tratio*j],1) Winc = np.sum(dW[:,np.arange(Tratio*(j-1)+1,Tratio*j)],1) #Xtemp = Xtemp + Dt*lamda*Xtemp + mu*Xtemp*Winc; Xtemp = Xtemp + self._drift(x=Xtemp) + self._sig(x=Xtemp) * Winc #Dt*lamda*Xtemp + mu*Xtemp*Winc; Xem[:,j] = Xtemp return Xem ''' R = 4; Dt = R*dt; L = N/R; % L EM steps of size Dt = R*dt Xem = zeros(1,L); % preallocate for efficiency Xtemp = Xzero; for j = 1:L Winc = sum(dW(R*(j-1)+1:R*j)); Xtemp = Xtemp + Dt*lambda*Xtemp + mu*Xtemp*Winc; Xem(j) = Xtemp; end ''' class ExactDiffusion(AffineDiffusion): '''Diffusion that has an exact integral representation this is currently mainly for geometric, log processes ''' def __init__(self): pass def exactprocess(self, xzero, nobs, ddt=1., nrepl=2): '''ddt : discrete delta t should be the same as an AR(1) not tested yet ''' t = np.linspace(ddt, nobs*ddt, nobs) #expnt = np.exp(-self.lambd * t) expddt = np.exp(-self.lambd * ddt) normrvs = np.random.normal(size=(nrepl,nobs)) #do I need lfilter here AR(1) ? if mean reverting lag-coeff<1 #lfilter does not handle 2d arrays, it does? inc = self._exactconst(expddt) + self._exactstd(expddt) * normrvs return signal.lfilter([1.], [1.,-expddt], inc) def exactdist(self, xzero, t): expnt = np.exp(-self.lambd * t) meant = xzero * expnt + self._exactconst(expnt) stdt = self._exactstd(expnt) return stats.norm(loc=meant, scale=stdt) class ArithmeticBrownian(AffineDiffusion): ''' :math:: dx_t &= \\mu dt + \\sigma dW_t ''' def __init__(self, xzero, mu, sigma): self.xzero = xzero self.mu = mu self.sigma = sigma def _drift(self, *args, **kwds): return self.mu def _sig(self, *args, **kwds): return self.sigma def exactprocess(self, nobs, xzero=None, ddt=1., nrepl=2): '''ddt : discrete delta t not tested yet ''' if xzero is None: xzero = self.xzero t = np.linspace(ddt, nobs*ddt, nobs) normrvs = np.random.normal(size=(nrepl,nobs)) inc = self._drift + self._sigma * np.sqrt(ddt) * normrvs #return signal.lfilter([1.], [1.,-1], inc) return xzero + np.cumsum(inc,1) def exactdist(self, xzero, t): expnt = np.exp(-self.lambd * t) meant = self._drift * t stdt = self._sigma * np.sqrt(t) return stats.norm(loc=meant, scale=stdt) class GeometricBrownian(AffineDiffusion): '''Geometric Brownian Motion :math:: dx_t &= \\mu x_t dt + \\sigma x_t dW_t $x_t $ stochastic process of Geometric Brownian motion, $\\mu $ is the drift, $\\sigma $ is the Volatility, $W$ is the Wiener process (Brownian motion). ''' def __init__(self, xzero, mu, sigma): self.xzero = xzero self.mu = mu self.sigma = sigma def _drift(self, *args, **kwds): x = kwds['x'] return self.mu * x def _sig(self, *args, **kwds): x = kwds['x'] return self.sigma * x class OUprocess(AffineDiffusion): '''Ornstein-Uhlenbeck :math:: dx_t&=\\lambda(\\mu - x_t)dt+\\sigma dW_t mean reverting process TODO: move exact higher up in class hierarchy ''' def __init__(self, xzero, mu, lambd, sigma): self.xzero = xzero self.lambd = lambd self.mu = mu self.sigma = sigma def _drift(self, *args, **kwds): x = kwds['x'] return self.lambd * (self.mu - x) def _sig(self, *args, **kwds): x = kwds['x'] return self.sigma * x def exact(self, xzero, t, normrvs): #TODO: aggregate over time for process with observations for all t # i.e. exact conditional distribution for discrete time increment # -> exactprocess #TODO: for single t, return stats.norm -> exactdist expnt = np.exp(-self.lambd * t) return (xzero * expnt + self.mu * (1-expnt) + self.sigma * np.sqrt((1-expnt*expnt)/2./self.lambd) * normrvs) def exactprocess(self, xzero, nobs, ddt=1., nrepl=2): '''ddt : discrete delta t should be the same as an AR(1) not tested yet # after writing this I saw the same use of lfilter in sitmo ''' t = np.linspace(ddt, nobs*ddt, nobs) expnt = np.exp(-self.lambd * t) expddt = np.exp(-self.lambd * ddt) normrvs = np.random.normal(size=(nrepl,nobs)) #do I need lfilter here AR(1) ? lfilter does not handle 2d arrays, it does? from scipy import signal #xzero * expnt inc = ( self.mu * (1-expddt) + self.sigma * np.sqrt((1-expddt*expddt)/2./self.lambd) * normrvs ) return signal.lfilter([1.], [1.,-expddt], inc) def exactdist(self, xzero, t): #TODO: aggregate over time for process with observations for all t #TODO: for single t, return stats.norm expnt = np.exp(-self.lambd * t) meant = xzero * expnt + self.mu * (1-expnt) stdt = self.sigma * np.sqrt((1-expnt*expnt)/2./self.lambd) from scipy import stats return stats.norm(loc=meant, scale=stdt) def fitls(self, data, dt): '''assumes data is 1d, univariate time series formula from sitmo ''' # brute force, no parameter estimation errors nobs = len(data)-1 exog = np.column_stack((np.ones(nobs), data[:-1])) parest, res, rank, sing = np.linalg.lstsq(exog, data[1:], rcond=-1) const, slope = parest errvar = res/(nobs-2.) lambd = -np.log(slope)/dt sigma = np.sqrt(-errvar * 2.*np.log(slope)/ (1-slope**2)/dt) mu = const / (1-slope) return mu, lambd, sigma class SchwartzOne(ExactDiffusion): '''the Schwartz type 1 stochastic process :math:: dx_t = \\kappa (\\mu - \\ln x_t) x_t dt + \\sigma x_tdW \\ The Schwartz type 1 process is a log of the Ornstein-Uhlenbeck stochastic process. ''' def __init__(self, xzero, mu, kappa, sigma): self.xzero = xzero self.mu = mu self.kappa = kappa self.lambd = kappa #alias until I fix exact self.sigma = sigma def _exactconst(self, expnt): return (1-expnt) * (self.mu - self.sigma**2 / 2. /self.kappa) def _exactstd(self, expnt): return self.sigma * np.sqrt((1-expnt*expnt)/2./self.kappa) def exactprocess(self, xzero, nobs, ddt=1., nrepl=2): '''uses exact solution for log of process ''' lnxzero = np.log(xzero) lnx = super(self.__class__, self).exactprocess(xzero, nobs, ddt=ddt, nrepl=nrepl) return np.exp(lnx) def exactdist(self, xzero, t): expnt = np.exp(-self.lambd * t) #TODO: check this is still wrong, just guessing meant = np.log(xzero) * expnt + self._exactconst(expnt) stdt = self._exactstd(expnt) return stats.lognorm(loc=meant, scale=stdt) def fitls(self, data, dt): '''assumes data is 1d, univariate time series formula from sitmo ''' # brute force, no parameter estimation errors nobs = len(data)-1 exog = np.column_stack((np.ones(nobs),np.log(data[:-1]))) parest, res, rank, sing = np.linalg.lstsq(exog, np.log(data[1:]), rcond=-1) const, slope = parest errvar = res/(nobs-2.) #check denominator estimate, of sigma too low kappa = -np.log(slope)/dt sigma = np.sqrt(errvar * kappa / (1-np.exp(-2*kappa*dt))) mu = const / (1-np.exp(-kappa*dt)) + sigma**2/2./kappa if np.shape(mu)== (1,): mu = mu[0] # TODO: how to remove scalar array ? if np.shape(sigma)== (1,): sigma = sigma[0] #mu, kappa are good, sigma too small return mu, kappa, sigma class BrownianBridge: def __init__(self): pass def simulate(self, x0, x1, nobs, nrepl=1, ddt=1., sigma=1.): nobs=nobs+1 dt = ddt*1./nobs t = np.linspace(dt, ddt-dt, nobs) t = np.linspace(dt, ddt, nobs) wm = [t/ddt, 1-t/ddt] #wmi = wm[1] #wm1 = x1*wm[0] wmi = 1-dt/(ddt-t) wm1 = x1*(dt/(ddt-t)) su = sigma* np.sqrt(t*(1-t)/ddt) s = sigma* np.sqrt(dt*(ddt-t-dt)/(ddt-t)) x = np.zeros((nrepl, nobs)) x[:,0] = x0 rvs = s*np.random.normal(size=(nrepl,nobs)) for i in range(1,nobs): x[:,i] = x[:,i-1]*wmi[i] + wm1[i] + rvs[:,i] return x, t, su class CompoundPoisson: '''nobs iid compound poisson distributions, not a process in time ''' def __init__(self, lambd, randfn=np.random.normal): if len(lambd) != len(randfn): raise ValueError('lambd and randfn need to have the same number of elements') self.nobj = len(lambd) self.randfn = randfn self.lambd = np.asarray(lambd) def simulate(self, nobs, nrepl=1): nobj = self.nobj x = np.zeros((nrepl, nobs, nobj)) N = np.random.poisson(self.lambd[None,None,:], size=(nrepl,nobs,nobj)) for io in range(nobj): randfnc = self.randfn[io] nc = N[:,:,io] #print nrepl,nobs,nc #xio = randfnc(size=(nrepl,nobs,np.max(nc))).cumsum(-1)[np.arange(nrepl)[:,None],np.arange(nobs),nc-1] rvs = randfnc(size=(nrepl,nobs,np.max(nc))) print('rvs.sum()', rvs.sum(), rvs.shape) xio = rvs.cumsum(-1)[np.arange(nrepl)[:,None],np.arange(nobs),nc-1] #print xio.shape x[:,:,io] = xio x[N==0] = 0 return x, N ''' randn('state',100) % set the state of randn T = 1; N = 500; dt = T/N; t = [dt:dt:1]; M = 1000; % M paths simultaneously dW = sqrt(dt)*randn(M,N); % increments W = cumsum(dW,2); % cumulative sum U = exp(repmat(t,[M 1]) + 0.5*W); Umean = mean(U); plot([0,t],[1,Umean],'b-'), hold on % plot mean over M paths plot([0,t],[ones(5,1),U(1:5,:)],'r--'), hold off % plot 5 individual paths xlabel('t','FontSize',16) ylabel('U(t)','FontSize',16,'Rotation',0,'HorizontalAlignment','right') legend('mean of 1000 paths','5 individual paths',2) averr = norm((Umean - exp(9*t/8)),'inf') % sample error ''' if __name__ == '__main__': doplot = 1 nrepl = 1000 examples = []#['all'] if 'all' in examples: w = Diffusion() # Wiener Process # ^^^^^^^^^^^^^^ ws = w.simulateW(1000, nrepl=nrepl) if doplot: plt.figure() tmp = plt.plot(ws[0].T) tmp = plt.plot(ws[0].mean(0), linewidth=2) plt.title('Standard Brownian Motion (Wiener Process)') func = lambda t, W: np.exp(t + 0.5*W) us = w.expectedsim(func, nobs=500, nrepl=nrepl) if doplot: plt.figure() tmp = plt.plot(us[0].T) tmp = plt.plot(us[1], linewidth=2) plt.title('Brownian Motion - exp') #plt.show() averr = np.linalg.norm(us[1] - np.exp(9*us[2]/8.), np.inf) print(averr) #print us[1][:10] #print np.exp(9.*us[2][:10]/8.) # Geometric Brownian # ^^^^^^^^^^^^^^^^^^ gb = GeometricBrownian(xzero=1., mu=0.01, sigma=0.5) gbs = gb.simEM(nobs=100, nrepl=100) if doplot: plt.figure() tmp = plt.plot(gbs.T) tmp = plt.plot(gbs.mean(0), linewidth=2) plt.title('Geometric Brownian') plt.figure() tmp = plt.plot(np.log(gbs).T) tmp = plt.plot(np.log(gbs.mean(0)), linewidth=2) plt.title('Geometric Brownian - log-transformed') ab = ArithmeticBrownian(xzero=1, mu=0.05, sigma=1) abs = ab.simEM(nobs=100, nrepl=100) if doplot: plt.figure() tmp = plt.plot(abs.T) tmp = plt.plot(abs.mean(0), linewidth=2) plt.title('Arithmetic Brownian') # Ornstein-Uhlenbeck # ^^^^^^^^^^^^^^^^^^ ou = OUprocess(xzero=2, mu=1, lambd=0.5, sigma=0.1) ous = ou.simEM() oue = ou.exact(1, 1, np.random.normal(size=(5,10))) ou.exact(0, np.linspace(0,10,10/0.1), 0) ou.exactprocess(0,10) print(ou.exactprocess(0,10, ddt=0.1,nrepl=10).mean(0)) #the following looks good, approaches mu oues = ou.exactprocess(0,100, ddt=0.1,nrepl=100) if doplot: plt.figure() tmp = plt.plot(oues.T) tmp = plt.plot(oues.mean(0), linewidth=2) plt.title('Ornstein-Uhlenbeck') # SchwartsOne # ^^^^^^^^^^^ so = SchwartzOne(xzero=0, mu=1, kappa=0.5, sigma=0.1) sos = so.exactprocess(0,50, ddt=0.1,nrepl=100) print(sos.mean(0)) print(np.log(sos.mean(0))) doplot = 1 if doplot: plt.figure() tmp = plt.plot(sos.T) tmp = plt.plot(sos.mean(0), linewidth=2) plt.title('Schwartz One') print(so.fitls(sos[0,:],dt=0.1)) sos2 = so.exactprocess(0,500, ddt=0.1,nrepl=5) print('true: mu=1, kappa=0.5, sigma=0.1') for i in range(5): print(so.fitls(sos2[i],dt=0.1)) # Brownian Bridge # ^^^^^^^^^^^^^^^ bb = BrownianBridge() #bbs = bb.sample(x0, x1, nobs, nrepl=1, ddt=1., sigma=1.) bbs, t, wm = bb.simulate(0, 0.5, 99, nrepl=500, ddt=1., sigma=0.1) if doplot: plt.figure() tmp = plt.plot(bbs.T) tmp = plt.plot(bbs.mean(0), linewidth=2) plt.title('Brownian Bridge') plt.figure() plt.plot(wm,'r', label='theoretical') plt.plot(bbs.std(0), label='simulated') plt.title('Brownian Bridge - Variance') plt.legend() # Compound Poisson # ^^^^^^^^^^^^^^^^ cp = CompoundPoisson([1,1], [np.random.normal,np.random.normal]) cps = cp.simulate(nobs=20000,nrepl=3) print(cps[0].sum(-1).sum(-1)) print(cps[0].sum()) print(cps[0].mean(-1).mean(-1)) print(cps[0].mean()) print(cps[1].size) print(cps[1].sum()) #Note Y = sum^{N} X is compound poisson of iid x, then #E(Y) = E(N)*E(X) eg. eq. (6.37) page 385 in http://ee.stanford.edu/~gray/sp.html #plt.show()