Эквивалент функции ручки в R - PullRequest
Новая
       27

Эквивалент функции ручки в R

1 голос
/ 05 марта 2020

Я перемещаю некоторый код из Matlab в R и сталкиваюсь с некоторыми трудностями в конкретном случае функции handle. Это мой код Matlab: 

function Application_ChFun
clear;close all;clc;warning('off');
StepsYr = 10;

%% --parameters-- %%
S0 = 1;   
r  = 0.0; 
t0 = 0;   
T2 = 5;

gamma = 0.5; 
kappa = 0.3; 
rho   = -0.6; 
vBar  = 0.05; 
v0    = 0.04; 

NoOfPaths = 5e4; 
NoOfSteps = StepsYr*T2;
%% --Define model-- %%
cf = @(u,T)ChFun(u, T, kappa,vBar,gamma,rho, v0, r);
Vc = @(t,x)MktFun(cf,t,x,log(S0));

% Define bump size
bump_T  = 1e-4;
bump_K  = @(T)1e-4;

% Define derivatives
dC_dT   = @(T,K) (Vc(T + bump_T,K) - Vc(T ,K)) /  bump_T;
dC_dK   = @(T,K) (Vc(T,K + bump_K(T)) - Vc(T,K - bump_K(T))) / (2 * bump_K(T));
d2C_dK2 = @(T,K) (Vc(T,K + bump_K(T)) + Vc(T,K-bump_K(T)) - 2*Vc(T,K)) / bump_K(T)^2;

t   = t0;
S = S0+zeros(NoOfPaths,1);

for i = 1:NoOfSteps

    if i==1
        t_adj = 1/NoOfSteps;
        t     = t_adj; 
    end       

    % AAA perfectly matches with the R equivalent, but AAB and AAC do not.
    AAA = dC_dT(t,S);
    AAB = dC_dK(t,S);
    AAC = d2C_dK2(t,S);
end 

function value = MktFun(cf,T,x,x0)
value = CM_Proxy(cf,T,x,x0);

function value = CM_Proxy(ChF,T,K,x0)
K(K<1e-5)=1e-5;

alpha   = 0.75; 
c       = 3e2;
N_CM    = 2^12;

eta    = c/N_CM;
b      = pi/eta;
u      = [0:N_CM-1]*eta;
lambda = 2*pi/(N_CM*eta);
i      = complex(0,1);

u_new = u-(alpha+1)*i; 
cf    = exp(i*u_new*x0).*ChF(u_new,T);
psi   = cf./(alpha^2+alpha-u.^2+i*(2*alpha+1)*u);

SimpsonW         = 3+(-1).^[1:N_CM]-[1,zeros(1,N_CM-1)];
SimpsonW(N_CM)   = 0;
SimpsonW(N_CM-1) = 1;
FFTFun           = exp(i*b*u).*psi.*SimpsonW;
payoff           = real(eta*fft(FFTFun)/3);
strike           = exp(-b:lambda:b-lambda);
payoff_specific  = spline(strike,payoff,K);
value = exp(-log(K)*alpha).*payoff_specific/pi;

function cf=ChFun(u, tau, kappa,vBar,gamma,rho, v0, r)
i     = complex(0,1);

D_1  = sqrt(((kappa -i*rho*gamma.*u).^2+(u.^2+i*u)*gamma^2));
g    = (kappa- i*rho*gamma*u-D_1)./(kappa-i*rho*gamma*u+D_1);    

C = (1/gamma^2)*(1-exp(-D_1*tau))./(1-g.*exp(-D_1*tau)).*(kappa-gamma*rho*i*u-D_1);
A = i*u*r*tau + kappa*vBar*tau/gamma^2 * (kappa-gamma*rho*i*u-D_1)-2*kappa*vBar/gamma^2*log((1-g.*exp(-D_1*tau))./(1-g));

cf = exp(A + C * v0);

, где MktFun - стандартная функция. Когда g=dC_dK(t,S) вызывается в первую очередь, оценивается bump_K(T), а затем Vc(T,K + bump_K(T)) и Vc(T,K-bump_K(T)).

В RI есть следующее: 

Application_ChFun <- function(){

  StepsYr = 10

  ## --parameters-- ##
  S0 = 1  
  r  = 0.0 
  t0 = 0   
  T2 = 5

  gamma = 0.5 
  kappa = 0.3 
  rho   = -0.6   
  vBar  = 0.05  
  v0    = 0.04  

  NoOfPaths = 5e4 
  NoOfSteps = StepsYr*T2

  ## --Define model-- ##
  cf <- function(u,T) ChFun(u,T,kappa,vBar,gamma,rho, v0, r)
  Vc <- function(t,x) MktFun(cf,t,x,log(S0))

  # Define bump size
  bump_T = 1e-4
  bump_K <- function(T) 1e-4

  # Define derivatives
  dC_dT   <- function(T,K) (Vc(T + bump_T,K) - Vc(T ,K)) /  bump_T
  dC_dK   <- function(T,K) (Vc(T,K + bump_K(T)) - Vc(T,K - bump_K(T))) / (2 * bump_K(T))
  d2C_dK2 <- function(T,K) (Vc(T,K + bump_K(T)) + Vc(T,K - bump_K(T)) - 2*Vc(T,K)) / bump_K(T)^2

  t = t0
  S = S0+rep(0,NoOfPaths)

  for (i in 1:NoOfSteps){

    t_real = t

    if (i==1){
      t_adj = 1/NoOfSteps;
      t     = t_adj 
    }

    # AAA perfectly matches with the R's equivalent. But AAB and AAC do not.
    AAA = dC_dT(t,S) 
    AAB = dC_dK(t,S)
    AAC = d2C_dK2(t,S)

  }
}


MktFun <- function(cf,T,x,x0){
  return(CM_Proxy(cf,T,x,x0))
}

CM_Proxy <- function(ChF,T,K,x0){

  K[K<1e-5] = 1e-5

  alpha = 0.75 
  c = 3e2
  N_CM = 2^12

  eta = c/N_CM
  b = pi/eta
  u = (0:(N_CM-1))*eta
  lambda = 2*pi/(N_CM*eta)
  i = complex(real = 0, imaginary = 1)

  u_new = u - (alpha+1)*i # European call option.
  cf = exp(i*u_new*x0)*ChF(u_new,T)
  psi = cf/(alpha^2+alpha-u^2+i*(2*alpha+1)*u)

  SimpsonW  = 3+(-1)^(1:N_CM)-c(1,rep(0,N_CM-1))
  SimpsonW[N_CM] = 0
  SimpsonW[N_CM-1] = 1
  FFTFun = exp(i*b*u)*psi*SimpsonW
  payoff = Re(eta*fft(FFTFun)/3)
  strike = exp(seq(-b,b-lambda,lambda))
  K = as.vector(K)
  payoff_specific = stinepack::stinterp(strike,payoff,K)
  value = exp(-log(K)*alpha)*payoff_specific$y/pi

  return(value)
}

ChFun <- function(u, tau, kappa,vBar,gamma,rho, v0, r){

  i = complex(real = 0, imaginary = 1)

  D_1 = sqrt(((kappa - i*rho*gamma*u)^2 + (u^2+i*u)*gamma^2))
  g = (kappa - i*rho*gamma*u - D_1) / (kappa - i*rho*gamma*u + D_1)   

  C = (1/gamma^2)*(1-exp(-D_1*tau))/(1-g*exp(-D_1*tau))*(kappa-gamma*rho*i*u-D_1)
  A = i*u*r*tau + kappa*vBar*tau/gamma^2 * (kappa-gamma*rho*i*u-D_1) +
                                                      -2*kappa*vBar/gamma^2*log((1-g*exp(-D_1*tau))/(1-g))

  cf = exp(A + C * v0)

  return(cf)
}

Проблема в том, что в в этом случае g=dC_dK(t,S) вызывает напрямую Vc вместо того, чтобы сначала вызвать bump_k(T). Может кто-нибудь предложить решение?

1 Ответ

0 голосов
/ 05 марта 2020

Порядок оценки функций не обязательно вывернут (как кажется, вы ожидаете) настолько, насколько это необходимо. R пытается сделать что-то лениво, поэтому, если вы включите дорогую операцию, на которую никогда не ссылаются, она не будет реализована.

Возьмите этот пример: 

f1 <- function(a) { message("f1"); a + 1; }
f2 <- function(b) { message("f2"); f1(b) + 2; }
f3 <- function(d) { message("f3"); f2(f1(d) + 3) / f2(f1(d) + 4); }
f3(2)
# f3
# f2
# f1
# f1
# f2
# f1
# f1
# [1] 0.9

Когда вызывается f3 следующие вызовы f2 Когда f2 вызывается впервые (с f1(d)+3), f2 вызывается с неоцененным аргументом. Как только f2 пытается использовать его b, только тогда он оценивается и вызывается f1.

Если я посмотрю на стек вызовов при первом вызове f1, то увидим: 

Browse[2]> where
where 1 at #1: f1(b)
where 2 at #1: f2(f1(d) + 3)
where 3 at #1: f3(2)

, показывающий порядок функций: сначала вызывается f3, затем f2, затем оттуда f1.

...