title: Equivalent of handle function in R
body: <p>I'm moving some code from Matlab to R and I'm facing some difficulties in a particular case of the handle function. This is my Matlab code:</p>
  
  <pre><code>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&lt;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);
  </code></pre>
  
  <p>where <code>MktFun</code> is a standard function. When <code>g=dC_dK(t,S)</code> is called first of all <code>bump_K(T)</code> is evaluated and then <code>Vc(T,K + bump_K(T))</code> and <code>Vc(T,K-bump_K(T))</code>.</p>
  
  <p>In R I have the following:</p>
  
  <pre><code>Application_ChFun &lt;- 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 &lt;- function(u,T) ChFun(u,T,kappa,vBar,gamma,rho, v0, r)
    Vc &lt;- function(t,x) MktFun(cf,t,x,log(S0))
  
    # Define bump size
    bump_T = 1e-4
    bump_K &lt;- function(T) 1e-4
  
    # Define derivatives
    dC_dT   &lt;- function(T,K) (Vc(T + bump_T,K) - Vc(T ,K)) /  bump_T
    dC_dK   &lt;- function(T,K) (Vc(T,K + bump_K(T)) - Vc(T,K - bump_K(T))) / (2 * bump_K(T))
    d2C_dK2 &lt;- 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 &lt;- function(cf,T,x,x0){
    return(CM_Proxy(cf,T,x,x0))
  }
  
  CM_Proxy &lt;- function(ChF,T,K,x0){
  
    K[K&lt;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 &lt;- 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)
  }
  </code></pre>
  
  <p>The problem is that in this case is that <code>g=dC_dK(t,S)</code> calls directly <code>Vc</code> instead of calling <code>bump_k(T)</code> first. Could someone please suggest a solution?</p>