期權
SVCJ (SVJJ) 達菲等。Matlab中的模型實現
我正在嘗試在 MATLAB中實現Duffie 等人提出的上述 SVCJ 模型。到目前為止沒有成功。它應該為普通（歐洲）電話定價。提供的參數，預期價格為：~6.0411。該函式拋出錯誤：
警告：已達到使用的最大間隔數限制。誤差的近似界限是 1.1e+03。積分可能不存在，或者可能難以在數值上逼近所要求的精度。
使用 quadv 而不是積分並沒有產生預期的結果，而是只返回了 NaN。
function price = SVCJDuffie2000Test(S0, K, V0, r, mus, muv, lambda, kappav, ... thetav, sigmav, sigmas, rho, rhoj, q, t, T) % this function should in theory calculate the analytical solution for % the SVCJ model by Duffie et al % % % % % T=0.25 % K = 100; % S0 = 85; % kappav=4; % lambda=4; % rhoj=-0.5; % rho=-0.5; % thetav=0.04; % r=0.05; % sigmas=0.06; % mu = 0.2380 % muv=0.02; % muJ = -0.04; % sigmav=0.1; % V0 = 0.25; % % % mus = log((1+mu)*(1-rhoj*muv))-0.5*(sigmas^2); % SVCJDuffie2000Test( S0, K, V0, r, mus, muv, lambda, kappav, thetav, sigmav, sigmas, rho, rhoj, 0,0,T ) % S0, V0, r, mu, muv, lambda, kappav, % thetav, sigmav, sigmas, rho, rhoj, y = log(S0); X0 = y; % nu, initial vola, nubar long run mean vola nu = V0; c = K; nubar = thetav; sigmav = sigmav; sigmacy = sigmas; mucy = mus; mucv = muv; rhoj = rhoj; rhobar = rho; kappav = kappav; zetabar = q; lambdac = lambda; % specific to SVCJ model lambdabar = lambda; r = r; % not needed for SVJJ, ie SVCJ sigmay = 0; lambday = 0; lambdav = 0; muv = 0; muy = 0; mubar = thetafunc(1,0) - 1 function retval = alphabar(tau,u) thetainter = lambdabar^-1*(lambdac*fc(u,tau)); retval = alpha0(tau, u) - lambdabar*tau*(1+mubar*u)+lambdabar*thetainter; end function retval = betabar(tau,u) a = u*(1-u); b = sigmav*rhobar*u-kappav; gamma = sqrt(b^2+a*sigmav^2); retval = a*(1-exp(-gamma*tau)) / (2*gamma-(gamma+b)*(1-exp(-gamma*tau))); end function retval = alpha0(tau, u) a = u*(1-u); b = sigmav*rhobar*u-kappav; gamma = sqrt(b^2+a*sigmav^2); retval = -r*tau+(r-zetabar)*u*tau-kappav*nubar*... (((gamma+b)/sigmav^2)*tau + ... (2/sigmav^2)*log(1-((gamma+b)/(2*gamma))*(1-exp(-gamma*tau)))); end function retval = fy(u,tau) retval = tau*exp(muy*u+1/2*sigmay^2*u^2); end function retval = fv(u,tau) a = u*(1-u); b = sigmav*rhobar*u-kappav; gamma = sqrt(b^2+a*sigmav^2); c_ = 1 - rhoj*mucv*u; retval = (gamma-b)/((gamma-b)+muv*a)*tau-(2*muv*a)/((gamma*c_)^2-(b*c_-muv*a)^2)... *log(1-((gamma+b)*c_-muv*a)/(2*gamma)*(1-exp(-gamma*tau))); end function retval = fc(u,tau) a = u*(1-u); b = sigmav*rhobar*u-kappav; c_ = 1 - rhoj*mucv*u; gamma = sqrt(b^2+a*sigmav^2); d_ = (gamma-b)/((gamma-b)*c_+mucv*a)*tau-(2*mucv*a)/((gamma*c_)^2-(b*c_-mucv*a)^2)... *log(1-((gamma+b)*c_-mucv*a)/(2*gamma*c_)*(1-exp(-gamma*tau))); retval = exp(mucy*u+sigmacy^2*u^2/2)*d_; end function retval = thetafunc(c1, c2) retval = lambdabar^-1*(lambdac*thetacfunc(c1,c2)); end function retval = thetayfunc(c) retval = exp(muy*c+1/2*sigmay^2*c^2); end function retval = thetavfunc(c) retval = 1/(1-muv*c); end function retval = thetacfunc(c1,c2) retval = exp(mucy*c1+1/2*sigmacy^2*c1^2) / ... (1-mucv*c2-rhoj*mucv*c1); end function retval = psi(y_, u, t, T) retval = exp(alphabar(T-t,u)+u*y_+ betabar(T-t,u)*nu); end % extended transform p.13 function retval = psichi(nu_, u, y, t, T) retval = exp(alphabar(T-t,u)+u*y+ betabar(T-t,u)*nu_); end function retval = Gab(a,b, y, X0_, T) integratefunc = @(nu_) imag(psi(complex(a,nu_*b),X0_, 0, T) * exp(complex(0,-nu_*y))) / nu_; retval = psi(a,X0_,0,T)/2 - 1/pi*integral(integratefunc, 0, Inf, 'ArrayValued',true); %retval = 1; end % depends on payoff function, see p. 6, 18 %aT = -r*T; %bT = 0; %d = 1; % could also be: aT = 0; bT = 1; d = -r*T; GbTplusdminusd = Gab(bT+d,-d,-log(c),X0,T) GbTminusd =Gab(bT,-d,-log(c),X0,T) price = exp(aT)*GbTplusdminusd-c*exp(aT)*GbTminusd; end論文中的功能：
劇院：p。23
alpha0：第 22 頁
字母欄：第 22 頁
測試欄：第 22 頁
fy，fv，fc：p。23
thetayfunc, thetacfunc, thetavfunc: p。22
狗：第 21 頁
加布：p。13
不是函式，而是取決於衍生類型：
GbTplusdminusd，GbTminusd，價格：p。18
lambdabar = lambda，根據：p。24
- 關於積分問題：你的被積函式是高度振蕩的，Matlab 的自適應正交不能很好地處理這樣的被積函式。一般來說，當 Matlab 的標準程序表現不佳時，我會推薦 Mathematica。在這種情況下，Levin 類型的方法會執行得更好。quadv 產生 NaN 值的原因是由於 0 處的奇異性。您必須從一個小的正值開始積分。
- 無論如何，您的程式碼中有一些拼寫錯誤。我試圖更正程式碼。集成現在執行良好，但我仍然沒有得到您期望的值（我已將 t=0、q=0、mu=0 用於未指定的參數；在 mus 公式中需要 mu）。我已在已更改的行中添加了 XXX 作為註釋。此外，我使所有操作向量都被重視以加快集成速度。
function price = SVCJDuffie2000Test(S0, K, V0, r, mus, muv, lambda, kappav, ... thetav, sigmav, sigmas, rho, rhoj, q, t, T) % this function should in theory calculate the analytical solution for % the SVCJ model by Duffie et al % % % % % T=0.25 % K = 100; % S0 = 85; % kappav=4; % lambda=4; % rhoj=-0.5; % rho=-0.5; % thetav=0.04; % r=0.05; % sigmas=0.06; % muv=0.02; % muJ = -0.04; % sigmav=0.1; % V0 = 0.25; % % % mus = log((1+mu)*(1-rhoj*muv))-0.5*(sigmas^2); % SVCJDuffie2000Test( S0, K, V0, r, mus, muv, lambda, kappav, thetav, sigmav, sigmas, rho, rhoj, 0,0,T ) % S0, V0, r, mu, muv, lambda, kappav, % thetav, sigmav, sigmas, rho, rhoj, y = log(S0); X0 = [y; V0]; % XXX two-dimensional process % nu, initial vola, nubar long run mean vola %nu = V0; % XXX c = K; nubar = thetav; %sigmav = sigmav; sigmacy = sigmas; mucy = mus; mucv = muv; %rhoj = rhoj; rhobar = rho; %kappav = kappav; zetabar = q; lambdac = lambda; % not needed for SVJJ, ie SVCJ sigmay = 0; lambday = 0; lambdav = 0; muv = 0; muy = 0; % specific to SVCJ model lambdabar = lambday + lambdav + lambdac; r = r; mubar = thetafunc(1,0) - 1; function retval = alphabar(tau,u) thetainter = lambdabar^-1*(lambday*fy(u,tau)+lambdav*fv(u,tau)+lambdac*fc(u,tau)); retval = alpha0(tau, u) - lambdabar*tau*(1+mubar*u)+lambdabar*thetainter; end function retval = betabar(tau,u) a = u.*(1-u); b = sigmav*rhobar*u-kappav; gamma = sqrt(b.^2+a*sigmav^2); retval = -a.*(1-exp(-gamma*tau)) ./ (2*gamma-(gamma+b).*(1-exp(-gamma*tau))); % XXX minus was missing end function retval = alpha0(tau, u) a = u.*(1-u); b = sigmav*rhobar*u-kappav; gamma = sqrt(b.^2+a*sigmav^2); retval = -r*tau+(r-zetabar)*u*tau-kappav*nubar*... (((gamma+b)/sigmav^2)*tau + ... (2/sigmav^2)*log(1-((gamma+b)./(2*gamma)).*(1-exp(-gamma*tau)))); end function retval = fy(u,tau) retval = tau*exp(muy*u+1/2*sigmay^2*u.^2); end function retval = fv(u,tau) a = u.*(1-u); b = sigmav*rhobar*u-kappav; gamma = sqrt(b.^2+a*sigmav^2); retval = (gamma-b)./((gamma-b)+muv*a)*tau-(2*muv*a)./(gamma.^2-(b-muv*a).^2)... .*log(1-((gamma+b)-muv*a)./(2*gamma).*(1-exp(-gamma*tau))); % XXX equation had a number of mistakes retval(a==0) = (gamma-b)./((gamma-b)+muv*a)*tau; % XXX take care of special case end function retval = fc(u,tau) a = u.*(1-u); b = sigmav*rhobar*u-kappav; c_ = 1 - rhoj*mucv*u; gamma = sqrt(b.^2+a*sigmav^2); d_ = (gamma-b)./((gamma-b).*c_+mucv*a)*tau-(2*mucv*a)./((gamma.*c_).^2-(b.*c_-mucv*a).^2)... .*log(1-((gamma+b).*c_-mucv*a)/(2*gamma.*c_).*(1-exp(-gamma*tau))); retval = exp(mucy*u+sigmacy^2*u.^2/2).*d_; retval(a==0) = (gamma-b)./((gamma-b).*c_+mucv*a)*tau; % XXX take care of special case end function retval = thetafunc(c1, c2) retval = lambdabar^-1*(lambday*thetayfunc(c1)+lambday*thetavfunc(c2)+lambdac*thetacfunc(c1,c2)); end function retval = thetayfunc(c) retval = exp(muy*c+1/2*sigmay^2*c^2); end function retval = thetavfunc(c) retval = 1/(1-muv*c); end function retval = thetacfunc(c1,c2) retval = exp(mucy*c1+1/2*sigmacy^2*c1^2) / ... (1-mucv*c2-rhoj*mucv*c1); end % equation 4.2, p. 21 function retval = psi(u, y, v, t, T) retval = exp(alphabar(T-t,u)+u*y+ betabar(T-t,u)*v); end function retval = psichi(u, X0, t, T) y = X0(1); v = X0(2); retval = psi(u, y, v, t, T); end function retval = Gab(a,b, y, X0_, T) integratefunc = @(nu_) imag(psichi(complex(a,nu_*b),X0_, 0, T) .* exp(complex(0,-nu_*y))) ./ nu_; % XXX retval = psichi(a,X0_,0,T)/2 - 1/pi*integral(integratefunc, 0, Inf, 'ArrayValued', true); %retval = 1; end % depends on payoff function, see p. 6, 18 aT = 0; % XXX values for European call option bT = 0; d = 1; % could also be: %aT = 0; %bT = 1; %d = -r*T; GbTplusdminusd = Gab(bT+d,-d,-log(c),X0,T) GbTminusd =Gab(bT,-d,-log(c),X0,T) price = exp(aT)*GbTplusdminusd-c*exp(aT)*GbTminusd; end