В этом примере сравниваются альтернативные реализации разделяемого многомерного геометрического броуновского процесса движения, который часто называют многомерной рыночной моделью. Он моделирует пути выборки портфеля индексов акций с использованием sde, sdeddo, sdeld, cev, и gbm объекты.
Моделируемая рыночная модель:
| (1) |
где:
λ - диагональная матрица ожидаемых результатов индекса.
D - диагональная матрица с Xt вдоль диагонали.
λ - диагональная матрица стандартных отклонений результатов индекса.
Создание sde объект с использованием sde для представления модели рынка акций.
Загрузить Data_GlobalIdx2 набор данных:
load Data_GlobalIdx2 prices = [Dataset.TSX Dataset.CAC Dataset.DAX ... Dataset.NIK Dataset.FTSE Dataset.SP];
Перевести дневные цены в возвраты:
returns = tick2ret(prices);
Вычислить статистику данных для ввода в методы моделирования:
nVariables = size(returns, 2); expReturn = mean(returns); sigma = std(returns); correlation = corrcoef(returns); t = 0; X = 100; X = X(ones(nVariables,1));
Создайте простые анонимные функции дрейфа и диффузии, доступные для (t, Xt):
F = @(t,X) diag(expReturn) * X; G = @(t,X) diag(X) * diag(sigma);
Используйте эти функции с sde для создания sde объект для представления рыночной модели в уравнении 1:
SDE = sde(F, G, 'Correlation', correlation, 'StartState', X)
SDE =
Class SDE: Stochastic Differential Equation
-------------------------------------------
Dimensions: State = 6, Brownian = 6
-------------------------------------------
StartTime: 0
StartState: 100 (6x1 double array)
Correlation: 6x6 double array
Drift: drift rate function F(t,X(t))
Diffusion: diffusion rate function G(t,X(t))
Simulation: simulation method/function simByEuler
sde объект требует дополнительной информации для определения размерности модели, поскольку функции переданыsde известны только по их интерфейсу (t, Xt). Другими словами,sde объект требует только двух входов: функции скорости дрейфа и функции скорости диффузии, которые доступны как путем прохождения времени выборки, так и соответствующего вектора состояния (t, Xt).
В этом случае этой информации недостаточно для однозначного определения размерности вектора состояния и броуновского движения. Размерность разрешается путем задания вектора начального состояния. StartState. Подсистема SDE назначила метод моделирования по умолчанию, simByEuler, в Simulation параметр.
Создание sdeddo объект с использованием sdeddo для представления рыночной модели в уравнении 1:
Загрузить Data_GlobalIdx2 набор данных:
load Data_GlobalIdx2 prices = [Dataset.TSX Dataset.CAC Dataset.DAX ... Dataset.NIK Dataset.FTSE Dataset.SP];
Перевести дневные цены в возвраты:
returns = tick2ret(prices);
Вычислить статистику данных для ввода в методы моделирования:
nVariables = size(returns, 2); expReturn = mean(returns); sigma = std(returns); correlation = corrcoef(returns);
Создать drift и diffusion объекты с использованием drift и diffusion:
F = drift(zeros(nVariables,1), diag(expReturn))
F =
Class DRIFT: Drift Rate Specification
-------------------------------------
Rate: drift rate function F(t,X(t))
A: 6x1 double array
B: 6x6 diagonal double array
G = diffusion(ones(nVariables,1), diag(sigma))
G =
Class DIFFUSION: Diffusion Rate Specification
---------------------------------------------
Rate: diffusion rate function G(t,X(t))
Alpha: 6x1 double array
Sigma: 6x6 diagonal double array
Пройти drift и diffusion объекты в sdeddo:
SDEDDO = sdeddo(F, G, 'Correlation', correlation, ... 'StartState', 100)
SDEDDO =
Class SDEDDO: SDE from Drift and Diffusion Objects
--------------------------------------------------
Dimensions: State = 6, Brownian = 6
--------------------------------------------------
StartTime: 0
StartState: 100 (6x1 double array)
Correlation: 6x6 double array
Drift: drift rate function F(t,X(t))
Diffusion: diffusion rate function G(t,X(t))
Simulation: simulation method/function simByEuler
A: 6x1 double array
B: 6x6 diagonal double array
Alpha: 6x1 double array
Sigma: 6x6 diagonal double array
sdeddo объект требует два входных объекта, которые предоставляют больше информации, чем две функции из шага 4 представления рыночных моделей с использованием объектов SDE. Таким образом, размерность легче поддается разрешению. Фактически, начальная цена каждого индекса является скаляром (StartState = 100). Это в отличие от sde , который требовал явного вектора состояния для однозначного определения размерности задачи.
Снова класс каждого объекта четко идентифицируется, и параметры отображаются как поля структуры. В частности, Rate параметр объектов дрейфа и диффузии идентифицируется как вызываемая функция времени и состояния, F (t, Xt) и G (t, Xt) соответственно. Дополнительные параметры,A, B, Alpha, и Sigma, представляют собой массивы соответствующей размерности, указывающие статические (не изменяющиеся во времени) параметры. Другими словами, A (t, Xt), B (t, Xt), Alpha (t, Xt) и Sigma (t, Xt) являются постоянными функциями времени и состояния.
Создать sdeld, cev, и gbm объекты для представления рыночной модели в уравнении 1.
Загрузить Data_GlobalIdx2 набор данных:
load Data_GlobalIdx2 prices = [Dataset.TSX Dataset.CAC Dataset.DAX ... Dataset.NIK Dataset.FTSE Dataset.SP];
Перевести дневные цены в возвраты:
returns = tick2ret(prices);
Вычислить статистику данных для ввода в методы моделирования:
nVariables = size(returns, 2); expReturn = mean(returns); sigma = std(returns); correlation = corrcoef(returns); t = 0; X = 100; X = X(ones(nVariables,1));
Создание sdeld объект с использованием sdeld:
SDELD = sdeld(zeros(nVariables,1), diag(expReturn), ... ones(nVariables,1), diag(sigma),'Correlation', ... correlation, 'StartState', X)
SDELD =
Class SDELD: SDE with Linear Drift
----------------------------------------
Dimensions: State = 6, Brownian = 6
----------------------------------------
StartTime: 0
StartState: 100 (6x1 double array)
Correlation: 6x6 double array
Drift: drift rate function F(t,X(t))
Diffusion: diffusion rate function G(t,X(t))
Simulation: simulation method/function simByEuler
A: 6x1 double array
B: 6x6 diagonal double array
Alpha: 6x1 double array
Sigma: 6x6 diagonal double array
Создать cev объект с использованием cev:
CEV = cev(diag(expReturn), ones(nVariables,1), ... diag(sigma), 'Correlation', correlation, ... 'StartState', X)
CEV =
Class CEV: Constant Elasticity of Variance
------------------------------------------
Dimensions: State = 6, Brownian = 6
------------------------------------------
StartTime: 0
StartState: 100 (6x1 double array)
Correlation: 6x6 double array
Drift: drift rate function F(t,X(t))
Diffusion: diffusion rate function G(t,X(t))
Simulation: simulation method/function simByEuler
Return: 6x6 diagonal double array
Alpha: 6x1 double array
Sigma: 6x6 diagonal double array
Создать gbm объект с использованием gbm:
GBM = gbm(diag(expReturn), diag(sigma), 'Correlation', ... correlation, 'StartState', X)
GBM =
Class GBM: Generalized Geometric Brownian Motion
------------------------------------------------
Dimensions: State = 6, Brownian = 6
------------------------------------------------
StartTime: 0
StartState: 100 (6x1 double array)
Correlation: 6x6 double array
Drift: drift rate function F(t,X(t))
Diffusion: diffusion rate function G(t,X(t))
Simulation: simulation method/function simByEuler
Return: 6x6 diagonal double array
Sigma: 6x6 diagonal double array
Обратите внимание на последовательность ограничений интерфейса:
Однако все три объекта представляют одну и ту же многомерную рыночную модель.
Также, cev и gbm объекты отображают базовый параметр B полученные из sdeld объект как Return. Это интуитивно понятное название, обычно связанное с моделями долевого участия.
Загрузить Data_GlobalIdx2 набор и использование данных sde чтобы указать модель SDE как в разделе Представление рыночных моделей с использованием объектов SDE.
load Data_GlobalIdx2 prices = [Dataset.TSX Dataset.CAC Dataset.DAX ... Dataset.NIK Dataset.FTSE Dataset.SP]; returns = tick2ret(prices); nVariables = size(returns,2); expReturn = mean(returns); sigma = std(returns); correlation = corrcoef(returns); t = 0; X = 100; X = X(ones(nVariables,1)); F = @(t,X) diag(expReturn)* X; G = @(t,X) diag(X) * diag(sigma); SDE = sde(F, G, 'Correlation', ... correlation, 'StartState', X);
Моделирование единого пути коррелированных цен индекса собственного капитала в течение одного календарного года (определенного как приблизительно 250 торговых дней) с использованием значения по умолчаниюsimulate способ:
nPeriods = 249; % # of simulated observations dt = 1; % time increment = 1 day rng(142857,'twister') [S,T] = simulate(SDE, nPeriods, 'DeltaTime', dt); whos S
Name Size Bytes Class Attributes S 250x6 12000 double
Выходной массив S равно 250 на 6 = (NPeriods + 1около-nVariablesоколо-1) массив с тем же начальным значением, 100, для всех индексов. Каждая строка S является наблюдением вектора состояния Xt в момент времени t.
Постройте график моделируемых путей.
plot(T, S), xlabel('Trading Day'), ylabel('Price') title('Single Path of Multi-Dimensional Market Model') legend({'Canada' 'France' 'Germany' 'Japan' 'UK' 'US'}, ... 'Location', 'Best')
Поскольку simByEuler является допустимым методом моделирования, вы можете вызвать его напрямую, переопределяя Simulation текущий метод или функция параметра (который в данном случае simByEuler). Следующие инструкции создают те же ценовые пути, что и при моделировании рынков акций с использованием метода моделирования по умолчанию:
Загрузить Data_GlobalIdx2 набор и использование данных sde чтобы указать модель SDE как в разделе Представление рыночных моделей с использованием объектов SDE.
load Data_GlobalIdx2 prices = [Dataset.TSX Dataset.CAC Dataset.DAX ... Dataset.NIK Dataset.FTSE Dataset.SP]; returns = tick2ret(prices); nVariables = size(returns,2); expReturn = mean(returns); sigma = std(returns); correlation = corrcoef(returns); t = 0; X = 100; X = X(ones(nVariables,1)); F = @(t,X) diag(expReturn)* X; G = @(t,X) diag(X) * diag(sigma); SDE = sde(F, G, 'Correlation', ... correlation, 'StartState', X);
Моделирование одного пути с помощью simByEuler.
nPeriods = 249; % # of simulated observations dt = 1; % time increment = 1 day rng(142857,'twister') [S,T] = simByEuler(SDE, nPeriods, 'DeltaTime', dt);
Смоделировать 10 испытаний с одинаковыми исходными условиями и изучить S:
rng(142857,'twister') [S,T] = simulate(SDE, nPeriods, 'DeltaTime', dt, 'nTrials', 10); whos S
Name Size Bytes Class Attributes S 250x6x10 120000 double
Теперь выходной массив S является NPeriods + 1около-nVariablesоколо-NTrials массив временных рядов.
Постройте первые контуры.
plot(T, S(:,:,1)), xlabel('Trading Day'), ylabel('Price') title('First Path of Multi-Dimensional Market Model') legend({'Canada' 'France' 'Germany' 'Japan' 'UK' 'US'},... 'Location', 'Best')
Первая реализация S идентичен контурам на графике.
Наконец, рассмотрите возможность моделирования с помощью GBM методы моделирования. Отделимый GBM модели имеют два специальных метода моделирования:
Перегруженный метод моделирования Эйлера, simByEuler, предназначен для оптимальной производительности
Функция, simBySolution, обеспечивает приблизительное решение лежащего в основе стохастического дифференциального уравнения, рассчитанного на точность
Загрузить Data_GlobalIdx2 набор и использование данных sde определение модели SDE как в представлении рыночных моделей с использованием объектов SDE, а модели GBM как в представлении рыночных моделей с использованием объектов SDELD, CEV и GBM.
load Data_GlobalIdx2 prices = [Dataset.TSX Dataset.CAC Dataset.DAX ... Dataset.NIK Dataset.FTSE Dataset.SP]; returns = tick2ret(prices); nVariables = size(returns,2); expReturn = mean(returns); sigma = std(returns); correlation = corrcoef(returns); t = 0; X = 100; X = X(ones(nVariables,1)); F = @(t,X) diag(expReturn)* X; G = @(t,X) diag(X) * diag(sigma); SDE = sde(F, G, 'Correlation', ... correlation, 'StartState', X); GBM = gbm(diag(expReturn),diag(sigma), 'Correlation', ... correlation, 'StartState', X);
Чтобы проиллюстрировать преимущество производительности перегруженного метода аппроксимации Эйлера, увеличьте число испытаний до 10000.
nPeriods = 249; % # of simulated observations dt = 1; % time increment = 1 day rng(142857,'twister') [X,T] = simulate(GBM, nPeriods, 'DeltaTime', dt, ... 'nTrials', 10000); whos X
Name Size Bytes Class Attributes X 250x6x10000 120000000 double
Продукция X является гораздо большим массивом временных рядов.
Используя этот размер выборки, изучите терминальное распределение TSX Composite в Канаде, чтобы качественно проверить логнормальный характер данных.
histogram(squeeze(X(end,1,:)), 30), xlabel('Price'), ylabel('Frequency') title('Histogram of Prices after One Year: Canada (TSX Composite)')
Смоделировать 10 испытаний решения и построить график первого испытания:
rng(142857,'twister') [S,T] = simulate(SDE, nPeriods, 'DeltaTime', dt, 'nTrials', 10); rng(142857,'twister') [X,T] = simBySolution(GBM, nPeriods,... 'DeltaTime', dt, 'nTrials', 10); subplot(2,1,1) plot(T, S(:,:,1)), xlabel('Trading Day'),ylabel('Price') title('1st Path of Multi-Dim Market Model:Euler Approximation') subplot(2,1,2) plot(T, X(:,:,1)), xlabel('Trading Day'),ylabel('Price') title('1st Path of Multi-Dim Market Model:Analytic Solution')
В этом примере все параметры являются константами, и simBySolution действительно пробует точное решение. Детали единого индекса для любого данного испытания показывают, что ценовые пути аппроксимации Эйлера и точного решения близки, но не идентичны.
Следующий график иллюстрирует разницу между двумя методами:
subplot(1,1,1) plot(T, S(:,1,1) - X(:,1,1), 'blue'), grid('on') xlabel('Trading Day'), ylabel('Price Difference') title('Euler Approx Minus Exact Solution:Canada(TSX Composite)')
simByEuler Аппроксимация Эйлера буквально оценивает стохастическое дифференциальное уравнение непосредственно из уравнения движения, для некоторого подходящего значения dt приращение времени. Эта простая аппроксимация страдает ошибкой дискретизации. Эта ошибка объясняется несоответствием между выбором dt приращение времени и то, что в теории является непрерывным параметром времени.
Дискретно-временная аппроксимация улучшается DeltaTime приближается к нулю. Метод Эйлера часто является наименее точным и наиболее общим доступным методом. Этот метод применяется ко всем моделям, поставляемым в комплект моделирования.
Напротив,simBySolution способ обеспечивает более точное описание базовой модели. Этот метод моделирует ценовые пути путем аппроксимации замкнутого решения разделяемых моделей. В частности, он применяет эйлеровский подход к преобразованному процессу, который в целом не является точным решением этого GBM модель. Это происходит потому, что распределения вероятности моделируемого и истинного векторов состояния идентичны только для кусочно-постоянных параметров.
Когда все параметры модели являются кусочно постоянными в течение каждого периода наблюдения, смоделированный процесс является точным для времени наблюдения, в которое дискретизируется вектор состояния. Поскольку в этом примере все параметры являются константами,simBySolution действительно пробует точное решение.
Пример использования simBySolution для оптимизации точности решений см. раздел Оптимизация точности: о точности решения и ошибке.
Этот пример иллюстрирует две методики, которые индуцируют зависимость между отдельными элементами вектора состояния. Он также иллюстрирует взаимодействие между Sigma и Correlation.
Первая методика генерирует коррелированные гауссовы вариации для формирования броуновского процесса движения с зависимыми компонентами. Затем эти компоненты взвешиваются по диагональной волатильности или матрице воздействия. Sigma.
Вторая методика генерирует независимые гауссовы вариации для формирования стандартного броуновского процесса движения, который затем взвешивается более низким коэффициентом Холески желаемой ковариационной матрицы. Хотя эти методы могут быть использованы на многих моделях, взаимосвязь между ними легче всего проиллюстрировать, работая с разделяемым GBM модель (см. Моделирование цен собственного капитала с использованием методов моделирования GBM). Моделируемая рыночная модель:
где λ - диагональная матрица.
Загрузите набор данных:
load Data_GlobalIdx2 prices = [Dataset.TSX Dataset.CAC Dataset.DAX ... Dataset.NIK Dataset.FTSE Dataset.SP];
Преобразуйте дневные цены в возвраты:
returns = tick2ret(prices);
Определить Sigma и Correlation с использованием первой методики:
Используя первый метод, укажите Sigma в качестве диагональной матрицы стандартных отклонений возврата основных средств:
expReturn = diag(mean(returns)); % expected return vector sigma = diag(std(returns)); % volatility of returns
Определить Correlation в качестве выборки корреляционной матрицы этих возвращений. В этом случае компоненты броуновского движения зависимы:
correlation = corrcoef(returns); GBM1 = gbm(expReturn,sigma,'Correlation',... correlation);
Определить Sigma и Correlation используя второй метод:
Используя второй метод, укажите Sigma в качестве нижнего коэффициента Холеского ковариационной матрицы возврата активов:
covariance = cov(returns); sigma = cholcov(covariance)';
Набор Correlation к единичной матрице:
GBM2 = gbm(expReturn,sigma);
Здесь, sigma фиксирует как корреляцию, так и величину неопределенности возврата активов. В отличие от первой техники, компоненты броуновского движения независимы. Кроме того, этот метод допускает назначение по умолчанию единичной матрицы Correlationи является более простым.
Смоделировать одно испытание из 1000 наблюдений (примерно четыре года ежедневных данных) с использованием обоих методов. По умолчанию все переменные состояния начинаются с 1:
rng(22814,'twister') [X1,T] = simByEuler(GBM1,1000); % correlated Brownian motion rng(22814,'twister') [X2,T] = simByEuler(GBM2,1000); % standard Brownian motion
При использовании одного и того же начального состояния случайного числа каждый метод генерирует идентичные пути цены актива:
subplot(2,1,1) plot(T, X1) title('Sample Paths from Correlated Brownian Motion') ylabel('Asset Price') subplot(2,1,2) plot(T, X2) title('Sample Paths from Standard Brownian Motion') xlabel('Trading Day') ylabel('Asset Price')
Как описано в разделе Создание объектов SDE, параметры объекта могут оцениваться так, как если бы они представляли собой функции MATLAB ®, доступные через общий интерфейс. Эта доступность обеспечивает впечатление динамического поведения независимо от того, действительно ли основные параметры изменяются во времени. Кроме того, поскольку параметры доступны через общий интерфейс, кажущиеся простыми линейные конструкции могут фактически представлять сложные нелинейные конструкции .
Например, рассмотрим одномерную модель геометрического броуновского движения (GBM) вида:
XtdWt
В этой модели доходность, λ (t), и волатильность, λ (t), являются динамическими параметрами только времени. Однако при создании gbm для представления базовой модели такое динамическое поведение должно быть доступно через общий интерфейс (t, Xt). Это отражает тот факт, что GBM модели (и другие) являются ограниченными параметризациями, производными от общих SDE класс.
Для удобства можно задать параметры ограниченных моделей, например: GBM модели, как традиционные массивы MATLAB соответствующего размера. В этом случае такие массивы представляют собой статический частный случай более общей динамической ситуации, доступной интерфейсу (t, Xt).
Кроме того, при вводе параметров в качестве функций конструкторы объектов могут проверить, что они возвращают массивы правильного размера, оценивая их в начальное время и состоянии. В противном случае конструкторы объектов не знают какой-либо конкретной функциональной формы.
Следующий пример иллюстрирует метод, который включает динамическое поведение путем отображения традиционного массива временных рядов MATLAB на вызываемую функцию с интерфейсом (t, Xt). Он также сравнивает функцию с идентичной моделью с постоянными параметрами.
Поскольку массивы временных рядов представляют динамическое поведение, которое должно захватываться функциями, доступными интерфейсу (t, Xt), необходимы утилиты для преобразования традиционных массивов временных рядов в вызываемые функции времени и состояния. В следующем примере показано, как это сделать с помощью функции преобразования ts2func (временной ряд для функционирования).
Загрузите данные. Загрузите ежедневный набор исторических данных, содержащий трехмесячные ставки Euribor и уровни индекса закрытия CAC 40 Франции, охватывающие временной интервал с 7 февраля 2001 года по 24 апреля 2006 года:
load Data_GlobalIdx2Моделирование нейтральных по риску путей выборки. Моделирование нейтральных для риска путей выборки индекса CAC 40 с использованием геометрического броуновского движения (GBM) модель:
startXtdWt
где r (t) представляет эволюцию безрисковой нормы прибыли.
Кроме того, предположим, что требуется ежегодная обработка соответствующей информации, полученной из ежедневных данных (ежегодная обработка данных необязательна, но полезна для сравнения с другими примерами), и что каждый календарный год включает 250 торговых дней:
dt = 1/250; returns = tick2ret(Dataset.CAC); sigma = std(returns)*sqrt(250); yields = Dataset.EB3M; yields = 360*log(1 + yields);
Сравните пути выборки из двух подходов к историческому моделированию, нейтральных по риску. Сравните результирующие пути выборки, полученные из двух подходов исторического моделирования, нейтральных к риску, где ежедневные доходности Euribor служат прокси для безрисковой нормы прибыли.
Первый подход определяет нейтральную по риску доходность как выборочную среднюю доходность Euribor и, следовательно, предполагает постоянную (нединамическую) безрисковую доходность:
nPeriods = length(yields); % Simulated observations rng(5713,'twister') obj = gbm(mean(yields),diag(sigma),'StartState',100)
obj =
Class GBM: Generalized Geometric Brownian Motion
------------------------------------------------
Dimensions: State = 1, Brownian = 1
------------------------------------------------
StartTime: 0
StartState: 100
Correlation: 1
Drift: drift rate function F(t,X(t))
Diffusion: diffusion rate function G(t,X(t))
Simulation: simulation method/function simByEuler
Return: 0.0278117
Sigma: 0.231906
[X1,T] = simulate(obj,nPeriods,'DeltaTime',dt);Напротив, второй подход определяет нейтральную к риску доходность как исторические временные ряды доходности Euribor. Поэтому она предполагает динамичную, но детерминированную норму прибыли; этот пример не иллюстрирует стохастические процентные ставки. Чтобы проиллюстрировать этот динамический эффект, используйте ts2func утилита:
r = ts2func(yields,'Times',(0:nPeriods - 1)');ts2func упаковывает указанный массив временных рядов в вызываемую функцию времени и состояния и синхронизирует его с необязательным вектором времени. Например:
r(0,100)
ans = 0.0470
вычисляет функцию при (t = 0, X t = 100) и возвращает первый наблюдаемый выход по Эврибору. Однако можно также вычислить результирующую функцию в любое промежуточное время t и состояние Xt:
r(7.5,200)
ans = 0.0472
Кроме того, следующая команда дает тот же результат при вызове только по времени:
r(7.5)
ans = 0.0472
Эквивалентность этих двух последних команд подчеркивает некоторые важные особенности.
При указании параметров в качестве функций они должны правильно вычисляться при прохождении скалярного, действительного времени выборки (t) и NVars-by-1 вектор состояния (Xt). Они также должны генерировать массив соответствующих размеров, который в первом случае является скалярной константой, а во втором - скалярной, кусочно-постоянной функцией только времени.
Использование времени (t) или состояния (Xt) не требуется. В текущем примере функция правильно вычисляется, когда прошло время, за которым следует состояние, тем самым удовлетворяя минимальным требованиям. Тот факт, что он также правильно вычисляется, когда прошло только время, просто указывает на то, что функция не требует вектора состояния Xt. Важно отметить, что он работает при прохождении (t, Xt).
Кроме того, ts2func функция выполняет интерполяцию кусочных констант с удержанием нулевого порядка (ZOH). Понятие кусочно-постоянных параметров широко распространено во всей архитектуре SDE и более подробно обсуждается в разделе Оптимизация точности: о точности решения и ошибке.
Выполните второе моделирование с использованием того же начального состояния случайных чисел. Выполните сравнение, выполнив второе моделирование с использованием того же начального состояния случайных чисел:
rng(5713,'twister') obj = gbm(r, diag(sigma),'StartState',100)
obj =
Class GBM: Generalized Geometric Brownian Motion
------------------------------------------------
Dimensions: State = 1, Brownian = 1
------------------------------------------------
StartTime: 0
StartState: 100
Correlation: 1
Drift: drift rate function F(t,X(t))
Diffusion: diffusion rate function G(t,X(t))
Simulation: simulation method/function simByEuler
Return: function ts2func/vector2Function
Sigma: 0.231906
X2 = simulate(obj,nPeriods,'DeltaTime',dt);Сравните два испытания моделирования. Постройте график ряда безрисковых эталонных показателей для сравнения двух испытаний моделирования:
subplot(2,1,1) plot(dates,100*yields) datetick('x') xlabel('Date') ylabel('Annualized Yield (%)') title('Risk Free Rate(3-Mo Euribor Continuously-Compounded)') subplot(2,1,2) plot(T,X1,'red',T,X2,'blue') xlabel('Time (Years)') ylabel('Index Level') title('Constant vs. Dynamic Rate of Return: CAC 40') legend({'Constant Interest Rates' 'Dynamic Interest Rates'},... 'Location', 'Best')
Пути близки, но не точны. Синяя линия на последнем сюжете использует все исторические данные Euribor, и иллюстрирует единственную пробную версию исторического моделирования.
Как описано в разделе Обеспечение положительных процентных ставок, все методы моделирования и интерполяции позволяют указать одну или несколько функций формы:
Xt)
для оценки в конце каждого времени выборки.
Соответствующий пример иллюстрирует простую общую функцию обработки на конец периода для обеспечения неотрицательных процентных ставок. В этом примере показана функция обработки, позволяющая полностью избежать выходных данных моделирования.
Рассмотрим ценообразование европейских опционов на акции с помощью моделирования Monte Carlo в рамках Black-Scholes-Merton. Предположим, что запас имеет следующие признаки:
Акции в настоящее время торгуются на уровне 100.
Акции не платят дивидендов.
Волатильность акций составляет 50% годовых.
Цена страйка опциона составляет 95.
Срок действия опции истекает через три месяца.
Безрисковая ставка постоянна и составляет 10% годовых.
Для решения этой задачи моделируют эволюцию базового запаса одномерной геометрической броуновской моделью движения (GBM) с постоянными параметрами:
5XtdWt
Кроме того, предположим, что цена акций моделируется ежедневно и каждый календарный месяц включает 21 торговый день:
strike = 95;
rate = 0.1;
sigma = 0.5;
dt = 1/252;
nPeriods = 63;
T = nPeriods*dt;
obj = gbm(rate,sigma,'StartState',100);Цель состоит в том, чтобы смоделировать независимые пути ежедневных котировок акций и рассчитать цену европейских опционов в качестве нейтрального по риску выборочного среднего показателя выплаты дисконтированного терминального опциона по истечении 63 дней с настоящего момента. В этом примере цены опционов вычисляются по двум подходам:
Моделирование Монте-Карло, которое явно запрашивает моделируемые пути заготовки в качестве выходных данных. Пути вывода затем используются для оценки опционов.
Функция обработки в конце периода, доступная по времени и состоянию, которая записывает цену терминального запаса для каждого пути образца. Эта функция обработки реализована как вложенная функция с доступом к совместно используемой информации. Дополнительные сведения см. в разделе Example_BlackScholes.m.
Перед моделированием вызовите файл примера для доступа к функции обработки в конце периода:
nTrials = 10000; % Number of independent trials (i.e., paths)
f = Example_BlackScholes(nPeriods,nTrials)f = struct with fields:
BlackScholes: @Example_BlackScholes/saveTerminalStockPrice
CallPrice: @Example_BlackScholes/getCallPrice
PutPrice: @Example_BlackScholes/getPutPrice
Моделирование 10000 независимых испытаний (пути выборки). Запросите смоделированные пути цены запаса в качестве выходного документа и определите функцию обработки на конец периода:
rng(88161,'twister') X = simBySolution(obj,nPeriods,'DeltaTime',dt,... 'nTrials',nTrials,'Processes',f.BlackScholes);
Вычислите цены опционов непосредственно по моделируемым путям цены акций. Поскольку это европейские опционы, игнорируйте все промежуточные цены на акции:
call = mean(exp(-rate*T)*max(squeeze(X(end,:,:)) - strike, 0))
call = 13.9342
put = mean(exp(-rate*T)*max(strike - squeeze(X(end,:,:)), 0))
put = 6.4166
Оцените опционы косвенно, вызвав вложенные функции:
f.CallPrice(strike,rate)
ans = 13.9342
f.PutPrice(strike,rate)
ans = 6.4166
Для справки, теоретические цены колла и пут, рассчитанные на основе опционных формул Блэка-Шоулза, являются 13.6953 и 6.3497соответственно.
Несмотря на то, что шаги 3 и 4 создают одинаковые цены опционов, последний подход работает непосредственно с ценами терминальных запасов каждого пути выборки. Поэтому он намного эффективнее памяти. В этом примере нет убедительной причины запрашивать выходные данные.
bates | bm | cev | cir | diffusion | drift | gbm | heston | hwv | interpolate | merton | sde | sdeddo | sdeld | sdemrd | simByEuler | simByQuadExp | simBySolution | simBySolution | simulate | ts2func