QuantAlgos

Продолжаем рассматривать алгоритмы построения улыбки волатильности. В этой статье будем находить "справедливые" цены опционов при помощи модели Хестона, которая относится к так называемым моделям стохастической волатильности. Хестон предложил использовать в качестве модели базового актива систему следующих уравнений:

$dS_t=r S_t dt+\sqrt{V_t}S_tdW_t^1$

$dV_t=k(\theta-V_t)dt+\sigma\sqrt{V_t}dW_t^2$

$Cov(dW_t^1dW_t^2)=\rho dt$,

где $S_t,V_t$- цена и волатильность базового актива соответственно, $dW_t^1,dW_t^2$- случайные броуновские процессы с корреляцией $\rho$.  $V_t$ - это квадратичный процесс с возвратом к среднему (mean reverting) со средним значением $\theta$ и интенсивностью k. $\sigma$ - среднеквадратичное отклонение волатильности, r - безрисковая ставка (для маржируемых равна 0, поэтому исключим сразу этот параметр для российского рынка).

Реальное статистическое распределение приращений цен базового актива  плохо соответствует Гауссовскому распределению, на основе которого была получена формула Блэка-Шоулза. Модель Хестона может описывать разные стат. распределения , например, коэффициент  $\rho$ может быть интерпретирован как корелляция между логарифмом приращения цены и волатильностью актива, что позволяет учитывать эффект "толстых хвостов" распределения. График плотности распределения приращения цены с разными значениями $\rho$ приведен в заглавии поста.

Цена европейского колл опциона для модели Хестона вычисляется по формуле:

$C(S_t,V_t,t,T)=S_tP_1-K\exp(-r(T-t))P_2$, где

$P_j(x,V_t,T,K)=1/2+\frac{1}{\pi}\int_0^\infty Re\Bigg(\frac{\exp(-i\phi\ln(K))f_j(x,V_t,T,\phi)}{i\phi}\Bigg)d\phi$

$x=\ln(S_t)$

$f_j(x,V_t,T,\phi)=\exp\{C(T-t,\phi)+D(T-t,\phi)V_t+i\phi x\}$

$C(T-t,\phi)=r\phi i r+\frac{a}{\sigma^2}[(b_j-\phi\sigma\phi i+d)\tau-2\ln(\frac{1-g\exp(dr)}{1-g})]$

$D(T-t,\phi)=\frac{b_j-\rho\sigma\phi i+d}{\sigma^2}\Bigg(\frac{1-\exp(dr)}{1-g\exp(dr)}\Bigg)$

$g=\frac{b_j-\rho\sigma\phi i+d}{b_j-\rho\sigma\phi i-d}$

$d=\sqrt{(\rho\sigma\phi i-b_j)^2-\sigma^2(2u_j\phi i-\phi^2)}$

для j=1,2, где

$u_1=1/2,u_2=-1/2,a=k\theta,b_1=k+\lambda-\rho\sigma, b_2=k+\lambda.$

Этот набор формул кажется сложным, однако решить их достаточно просто с помощью программы на C#, которая будет приведена ниже. Сложность составляет только вычисление интеграла с верхним бесконечным пределом в формуле для $P_j$, который находится с помощью числового метода Гаусса-Лагендре в той же программе. Также, для упрощения, можно сократить число параметров, убрав из них меру риска $\lambda$ , применив риск-нейтральный подход. В этом случае:

$a=k^*\theta^*,b_1=k^*-\rho\sigma,b_2=k^*$

Функция расчета формулы Хестона на языке C#:

//a-нижний предел интеграла (равен 0)
//b - верхний предел интеграла. Выбирается значение от 100 до 200, в зависимости от нужной точности
//delta - вычисляется грек дельта, который равен значению Р1 в формуле Хестона
double HestonCallGaussLegendre(double S,double K,double T,double r,double kappa,double theta,
							   double sigma,double lambda,double v0,double rho,int trap,
							   double a, double b,ref double delta)
        {
	        // Числовое интегрирование
	        double[] int1=new double[32];
            double[] int2 = new double[32];
	        double y;
	        for (int k=0; k< =31; k++) {
                y = (a + b) / 2.0 + (b - a) / 2.0 * X[k];
		        int1[k] = W[k]*HestonCF(y,S,K,r,T,rho,kappa,theta,lambda,sigma,v0,1,trap);
		        int2[k] = W[k]*HestonCF(y,S,K,r,T,rho,kappa,theta,lambda,sigma,v0,2,trap);
	        }

	        // Векторы для интегральной суммы
	        double I1 = VectorSum(int1);
	        double I2 = VectorSum(int2);

	        // Определение Р1 и Р2
	        double P1 = 0.5 + 1.0/Math.PI*I1*(b-a)/2;
	        double P2 = 0.5 + 1.0/Math.PI*I2*(b-a)/2;
            delta = P1;

	        // Цена колл опциона
	        return S*P1 - K*P2;
        }
// Функция суммирования элементов вектора
        double VectorSum(double[] A)
        {
            double sum = 0;
            double n = A.Length;
            for (int i = 0; i <= n - 1; i++)
                sum += A[i];
            return sum;
        }

private double HestonCF(Complex phi, double Spot, double Strike,
                double Rate, double T, double Rho, double Kappa,
                double Theta, double Lambda, double Sigma, double V,
                int Pnum, int trap)
        {
	        Complex		 S=new Complex(Spot  , 0.0);	// Цена базового актива
	        Complex		 K=new Complex(Strike, 0.0);	// Страйк
	        Complex	     r=new Complex(Rate  , 0.0);	// Безрисковая ставка (для марж. опционов =0)
	        Complex	   tau=new Complex(T     , 0.0);	// Период до экспирации в долях года
	        Complex      i=new Complex(0.0   , 1.0);	// Мнимая переменная
	        Complex    rho=new Complex(Rho   , 0.0);	// Параметр Хестона: Корреляция
            Complex kappa = new Complex(Kappa, 0.0);	// Параметр Хестона: Скорость возврата
            Complex theta = new Complex(Theta, 0.0);	// Параметр Хестона: уровень возвратности
            Complex lambda = new Complex(Lambda, 0.0);	// Параметр Хестона: мера риска (равна 0 для риск-нейтрального подхода)
            Complex sigma = new Complex(Sigma, 0.0);	// Параметр Хестона: Среднеквадратичное волатильности
            Complex v0 = new Complex(V, 0.0);	        // Параметр Хестона: Текущая волатильность
	        Complex	   two=new Complex(2.0   , 0.0);	// число 2 в комплексной форме
            Complex one = new Complex(1.0, 0.0);	    // число 1 в комплексной форме
	        Complex y, a, u, b, sigma2, d, g, G, C, D, c, f;
  	        y = rho*sigma*phi*i;
  	        a = kappa*theta;
	        if (Pnum==1) 
            {
		        // Первая характеристическая функция
		        u = 0.5;
	  	        b = kappa + lambda - rho*sigma;
	        }
	        else 
            {
		        // Вторая характеристическая функция
		        u = -0.5;
		        b = kappa + lambda;
	        }
  	        sigma2 = Complex.Pow(sigma,2);
  	        d = Complex.Sqrt((y-b)*(y-b) - sigma2*(two*u*phi*i - phi*phi));
  	        g = (b - y + d)/(b - y - d);
	        if (trap==1) 
            {
                // Версия модели "Little Heston Trap"
		        c = one/g;
		        G = (one - c*Complex.Exp(-d*tau))/(one - c);
		        C = r*i*phi*tau + a/sigma2*((b - rho*sigma*i*phi - d)*tau - two*Complex.Log(G));
                D = (b - rho * sigma * i * phi - d) / sigma2 * ((one - Complex.Exp(-d * tau)) / (one - c * Complex.Exp(-d * tau)));
	        }
	        else
            {
		        // Оригинальный вариант Хестона
                G = (one - g * Complex.Exp(d * tau)) / (one - g);
		        C = r*i*phi*tau + a/sigma2*((b - rho*sigma*i*phi + d)*tau - two*Complex.Log(G));
                D = (b - rho * sigma * i * phi + d) / sigma2 * ((one - Complex.Exp(d * tau)) / (one - g * Complex.Exp(d * tau)));
	        }
	        f = Complex.Exp(C + D*v0 + i*phi*Complex.Log(S));	

	        // Вычисление реальной части подинтегрального выражения
	        return (Complex.Exp(-i*phi*Complex.Log(K))*f/i/phi).Real;
        }

Следующий шаг - определение пяти параметров $k,\theta,\sigma,\rho,V$(V - текущая волатильность). Для этого нужно откалибровать модель по наблюдаемым рыночным ценам опционов. Применяем стандартный метод - берем выборку цен для опционов разных страйков за определенный период времени (вместе со сроками до экспирации), при этом рыночной ценой опциона считаем среднюю цену между бидом и аском ($C^M(K_i,T_i)=(bid_i+ask_i)/2$, и минимизируем следующее выражение, применяя нелинейный метод наименьших квадратов (МНК):

$\min S(\Omega)=\min_{\Omega}\sum_{i=1}^{N}w_i[C_i^{\Omega}(K_i,T_i)-C_i^M(K_i,T_i)]^2\leq\sum_{i=0}^Nw_i[bid_i-ask_i]^2$

где $\Omega$ - вектор параметров, $w_i$ - задаваемые веса (их выбор обсудим позже), N - размер выборки. Выражение в правой части означает,что полученные значения должны попадать в промежуток между бидом и аском наблюдаемых рыночных цен. Это ограничение, равно как и условие $2k\theta>\sigma^2$ (волатильность не может падать до 0) позволяет сузить диапазон решений, полученных с помощью МНК. Таких решений может быть несколько из-за того, что МНК, попадая в локальный минимум выражения $S(\Omega)$, останавливается и выдает не оптимальные значения. Таким образом, нахождение оптимальных параметров модели Хестона является нетривиальной задачей, и применяются следующие способы ее решения:

• сделать множество вычислений с помощью МНК, задавая различные значения начальных параметров, а затем выбрать минимальное из всех полученых $S(\Omega)$, получив таким образом соответствующие ему параметры модели;
• применить алгоритмы поиска глобального минимума, такие, как Differential Evolution, ASA и т.п. Недостаток таких алгоритмов в значительном времени, требуемом для нахождения параметров.

Веса можно задать в соответствии с формулой : $w_i=1/|bid_i-ask_i|$. Это интуитивный выбор, основанный на том, что, чем шире спред, тем больше свобода выбора в значении цены опциона. Для российского рынка лучшая аппроксимация получалась у меня при выборе одинакового значения весов, равного 1, но я не брал в рассмотрение слишком дальние страйки.

Получив параметры модели Хестона, мы сможем вычислить цены опционов $C(K_i,T_i)$ для любого страйка и периода до экспирации. Для наглядности мы сможем построить улыбку волатильности по значениям подразумеваемой волатильности из формулы Блэка-Шоулза, подставив в нее хестоновские цены опционов - см. график в начале поста.

Модель Хестона отражает реальное статистическое распределение приращений цены базового актива значительно лучше, чем это делает модель Блэка-Шоулза, в чем вы сможете убедиться, сравнивая реальные рыночные цены опционов с полученными по этой модели. Однако у нее есть один существенный недостаток, который проявляется в том, что, если до экспирации остается небольшой срок (около недели для российского рынка) цены крайних страйков модель определяет неверно, в терминах подразумеваемой волатильности - хвосты улыбки начинают расходиться:

Чтобы устранить этот недостаток мы должны перейти к применению модифицированной модели Хестона - модели Бэйтса, являющейся одной из лучших аппроксимаций, позволяющих с макимальной точностью находить "справедливые" цены опционов. Ее мы рассмотрим в следующей части цикла статей про улыбку волатильности.