Brownian interpolation of stochastic differential equations
[XT, T] = interpolate(MDL, T, Paths)
[XT, T] = interpolate(MDL, Paths, 'Name1', Value1,
'Name2', Value2, ...)
All classes in SDE Class Hierarchy.
This method performs a Brownian interpolation into a userspecified time series array, based on a piecewiseconstant Euler sampling approach.
Consider a vectorvalued SDE of the form:
$$d{X}_{t}=F(t,{X}_{t})dt+G(t,{X}_{t})d{W}_{t}$$
where:
X is an NVARSby1
state
vector.
F is an NVARSby1
driftrate
vectorvalued function.
G is an NVARSbyNBROWNS diffusionrate matrixvalued function.
W is an NBROWNSby1
Brownian
motion vector.
Given a userspecified time series array associated with this equation, this method performs a Brownian (stochastic) interpolation by sampling from a conditional Gaussian distribution. This sampling technique is sometimes called a Brownian bridge.
Note:
Unlike simulation methods, the 
MDL  Stochastic differential equation model. 
T  NTIMES element vector of interpolation times.
The length of this vector determines the number of rows in the interpolated
output time series XT . 
Paths  NPERIODS byNVARS byNTRIALS time
series array of sample paths of correlated state variables. For a
given trial, each row of this array is the transpose of the state
vector X_{t} at time t. Paths is
the initial time series array into which interpolate performs
the Brownian interpolation. 
Specify optional input arguments as variablelength lists of
matching parameter name/value pairs: 'Name1'
, Value1
, 'Name2'
, Value2
,
... and so on. The following rules apply when specifying parametername
pairs:
Specify the parameter name as a character vector, followed by its corresponding parameter value.
You can specify parameter name/value pairs in any order.
Parameter names are case insensitive.
You can specify unambiguous partial character vector matches.
Valid parameter names are:
Times  Vector of monotonically increasing observation times associated
with the time series input Paths . If you do not
specify a value for this parameter, Times is a
zerobased, unitincrement column vector of length NPERIODS .

Refine  Scalar logical flag that indicates whether interpolate uses
the interpolation times you request (see T ) to
refine the interpolation as new information becomes available. If
you do not specify a value for this argument or set it to FALSE (the
default value), interpolate bases the interpolation
only on the state information specified in Paths .
If you set Refine to TRUE , interpolate inserts
all new interpolated states into the existing Paths array
as they become available. This refines the interpolation grid available
to subsequent interpolation times during the current trial. 
Processes  Function or cell array of functions that indicates a sequence of background processes or state vector adjustments of the form $${X}_{t}=P(t,{X}_{t})$$ interpolate method
runs processing functions at each interpolation time. They must accept
the current interpolation time t, and the current
state vector X_{t}, and return
a state vector that may be an adjustment to the input state. If
you specify more than one processing function, If
you do not specify a processing function, 
XT  NTIMES byNVARS byNTRIALS time
series array of interpolated state variables. For a given trial, each
row of this array is the transpose of the interpolated state vector X_{t} at
time t. XT is the interpolated
time series formed by interpolating into the input Paths time
series array. 
T  NTIMES by1 column vector of interpolation
times associated with the output time series XT .
If the input interpolation time vector T contains
no missing observations (NaN s), this output is
the same time vector as T , but with the NaN s
removed. This reduces the length of T and the number
of rows of XT . 
Many applications require knowledge of the state vector at intermediate sample times that are initially unavailable. One way to approximate these intermediate states is to perform a deterministic interpolation. However, deterministic interpolation techniques fail to capture the correct probability distribution at these intermediate times. Brownian (or stochastic) interpolation captures the correct joint distribution by sampling from a conditional Gaussian distribution. This sampling technique is sometimes referred to as a Brownian Bridge.
The default stochastic interpolation technique is designed to interpolate into an existing time series and ignore new interpolated states as additional information becomes available. This technique is the usual notion of interpolation, which is called Interpolation without refinement.
Alternatively, the interpolation technique may insert new interpolated states into the existing time series upon which subsequent interpolation is based, by that means refining information available at subsequent interpolation times. This technique is called interpolation with refinement.
Interpolation without refinement is a more traditional technique, and is most useful when the input series is closely spaced in time. In this situation, interpolation without refinement is a good technique for inferring data in the presence of missing information, but is inappropriate for extrapolation. Interpolation with refinement is more suitable when the input series is widely spaced in time, and is useful for extrapolation.
The stochastic interpolation method is available to any model.
It is best illustrated, however, by way of a constantparameter Brownian
motion process. Consider a correlated, bivariate Brownian motion (BM
)
model of the form:
$$\begin{array}{l}d{X}_{1t}=0.3dt+0.2d{W}_{1t}0.1d{W}_{2t}\\ d{X}_{2t}=0.4dt+0.1d{W}_{1t}0.2d{W}_{2t}\\ E[d{W}_{1t}d{W}_{2t}]=\rho dt=0.5dt\end{array}$$
Create a bm
object to represent
the bivariate model:
mu = [0.3; 0.4];
sigma = [0.2 0.1; 0.1 0.2];
rho = [1 0.5; 0.5 1];
obj = bm(mu,sigma,'Correlation',rho);
Assuming that the drift (Mu
) and
diffusion (Sigma
) parameters are annualized, simulate
a single Monte Carlo trial of daily observations for one calendar
year (250 trading days):
rng default % make output reproducible dt = 1/250; % 1 trading day = 1/250 years [X,T] = simulate(obj,250,'DeltaTime',dt);
It is helpful to examine a small interval in detail.
Interpolate into the simulated time series with a Brownian bridge:
t = ((T(1) + dt/2):(dt/2):(T(end)  dt/2));
x = interpolate(obj,t,X,'Times',T);
Plot both the simulated and interpolated values:
plot(T,X(:,1),'.r',T,X(:,2),'.b') grid on; hold on; plot(t,x(:,1),'or',t,x(:,2),'ob') hold off; xlabel('Time (Years)') ylabel('State') title('BiVariate Brownian Motion: \rho = 0.5') axis([0.4999 0.6001 0.25 0.4])
In this plot:
The solid red and blue dots indicate the simulated states of the bivariate model.
The straight lines that connect the solid dots indicate intermediate states that would be obtained from a deterministic linear interpolation.
Open circles indicate interpolated states.
Open circles associated with every other interpolated state encircle solid dots associated with the corresponding simulated state. However, interpolated states at the midpoint of each time increment typically deviate from the straight line connecting each solid dot.
You can gain additional insight into the behavior of stochastic interpolation by regarding a Brownian bridge as a Monte Carlo simulation of a conditional Gaussian distribution.
This example examines the behavior of a Brownian bridge over a single time increment.
Divide a single time increment of length dt
into
10 subintervals:
mu = [0.3; 0.4]; sigma = [0.2 0.1; 0.1 0.2]; rho = [1 0.5; 0.5 1]; obj = bm(mu,sigma,'Correlation',rho); rng default; % make output reproducible dt = 1/250; % 1 trading day = 1/250 years [X,T] = simulate(obj,250,'DeltaTime',dt); n = 125; % index of simulated state near middle times = (T(n):(dt/10):T(n + 1)); nTrials = 25000; % # of Trials at each time
In each subinterval, take 25000 independent draws from a Gaussian distribution, conditioned on the simulated states to the left, and right:
average = zeros(length(times),1); variance = zeros(length(times),1); for i = 1:length(times) t = times(i); x = interpolate(obj,t(ones(nTrials,1)),... X,'Times',T); average(i) = mean(x(:,1)); variance(i) = var(x(:,1)); end
Plot the sample mean and variance of each state variable:
Note: The following graph plots the sample statistics of the first state variable only, but similar results hold for any state variable. 
subplot(2,1,1); hold on; grid on; plot([T(n) T(n + 1)],[X(n,1) X(n + 1,1)],'.b') plot(times, average, 'or') hold off; title('Brownian Bridge without Refinement: Sample Mean') ylabel('Mean') limits = axis; axis([T(n) T(n + 1) limits(3:4)]); subplot(2,1,2) hold on; grid on; plot(T(n),0,'.b',T(n + 1),0,'.b') plot(times, variance, '.r') hold('off'); title('Brownian Bridge without Refinement: Sample Variance') xlabel('Time (Years)') ylabel('Variance') limits = axis; axis([T(n) T(n + 1) limits(3:4)]);
The Brownian interpolation within the chosen interval, dt, illustrates the following:
The conditional mean of each state variable lies on a straightline segment between the original simulated states at each endpoint.
The conditional variance of each state variable is a quadratic function. This function attains its maximum midway between the interval endpoints, and is zero at each endpoint.
The maximum variance, although dependent upon the
actual model diffusionrate function G(t,X), is
the variance of the sum of NBROWNS
correlated Gaussian
variates scaled by the factor dt/4.
The previous plot highlights interpolation without refinement, in that none of the interpolated states take into account new information as it becomes available. If you had performed interpolation with refinement, new interpolated states would have been inserted into the time series and made available to subsequent interpolations on a trialbytrial basis. In this case, all random draws for any given interpolation time would be identical. Also, the plot of the sample mean would exhibit greater variability, but would still cluster around the straightline segment between the original simulated states at each endpoint. The plot of the sample variance, however, would be zero for all interpolation times, exhibiting no variability.
AitSahalia, Y. "Testing ContinuousTime Models of the Spot Interest Rate." The Review of Financial Studies, Spring 1996, Vol. 9, No. 2, pp. 385–426.
AitSahalia, Y. "Transition Densities for Interest Rate and Other Nonlinear Diffusions." The Journal of Finance, Vol. 54, No. 4, August 1999.
Glasserman, P. Monte Carlo Methods in Financial Engineering. New York, SpringerVerlag, 2004.
Hull, J. C. Options, Futures, and Other Derivatives, 5th ed. Englewood Cliffs, NJ: Prentice Hall, 2002.
Johnson, N. L., S. Kotz, and N. Balakrishnan. Continuous Univariate Distributions. Vol. 2, 2nd ed. New York, John Wiley & Sons, 1995.
Shreve, S. E. Stochastic Calculus for Finance II: ContinuousTime Models. New York: SpringerVerlag, 2004.