Radial Basis Function

[treeleaf]

Hi,

For my engineering school, we have a project on radial basis functions to do.
The teacher gave us the article: Option Pricing with Radial Basis Functions: A Tutorial: Wilmott Magazine Article by Alonso Peña.
The project is to price options using radial basis functions.

I code in C#. With XLW, I can export the C# functions in Excel.
I use dnAnalytics for matrix operations.

For the moment, I try to price an European call.

The problem I'm facing is that I get wrong results.

I've taken these parameters:

• N= 30
• T =1
• M =30
• a =0
• b =3,912023005
• c =0,53958938
• E =15
• r =0,05
• sigma= 0,3
• theta =0,5
• S =15

And I got the price of the call= -1,524E+189.

The problem is that all the students got bad results. The teacher took the time to see my code and said it was correct.


Here's the C# code:

        public static double g(double x, double E, string CallPutFlag)
        {

            // Option Payoff
            // IN: x, E, CallPutFlag
            // OUT: payoff en fonction de CallPutFlag
 
            if (CallPutFlag == "C") return Math.Max(Math.Exp(x) - E, 0);
            else return Math.Max(E - Math.Exp(x), 0);
        }

        public static double phi(double x, double c)
        {
            return Math.Sqrt(Math.Pow(x, 2) + Math.Pow(c, 2));
        }

        public static double phi_prime(double x, double c)
        {
            return x / phi(x, c);
        }

        public static double phi_seconde(double x, double c)
        {
            return 1 / phi(x, c)-Math.Pow(x,2)/Math.Pow(phi(x,c),3);
        }

        public static double alpha(double t)
        {
            return 0;
        }

        public static double beta(double t, double S, double r, double T, double E)
        {
            return S-E*Math.Exp(-r*(T-t));
        }
 
        [ExcelExport("RBF function")]
        public static double RBF(
            [Parameter("N")] int N,
            [Parameter("T")] double T,
            [Parameter("M")] int M,
            [Parameter("a")] double a,
            [Parameter("b")] double b,
            [Parameter("c")] double c,
            [Parameter("E")] double E,
            [Parameter("C ou P")] string CallPutFlag,
            [Parameter("r")] double r,
            [Parameter("sigma")] double sigma,
            [Parameter("theta")] double theta,
            [Parameter("S")] double S)
        {
            Vector ksi = new DenseVector(N);
            Matrix Phi = new DenseMatrix(N, N);
            Matrix Phi1 = new DenseMatrix(N, N);
            Matrix Phi2 = new DenseMatrix(N, N);
            Vector lambda = new DenseVector(N);
            Vector gi = new DenseVector(N);
            double deltaT = T / M, t = 0, Prix = 0, x=Math.Log(S);
            int i, j;
   
            // Création du vecteur ksi
            for (i = 0; i < N; i++)
            {
                ksi[i] = a + i * ((b - a) / (N - 1));
                gi[i] = g(ksi[i], E, CallPutFlag);
            }

            // Création de la matrice Phi, Phi' et Phi''
            for (i = 0; i < N; i++)
            {
                for (j = 0; j < N; j++)
                {
                    Phi[i, j] = phi(ksi[i] - ksi[j], c);
                    Phi1[i, j] = phi_prime(ksi[i] - ksi[j], c);
                    Phi2[i, j] = phi_seconde(ksi[i] - ksi[j], c);
                }
            }

            Matrix Temp = new DenseMatrix(N, N);
            // Résolution de Phi*lambda=g(ksi)

            Temp = Phi.Inverse();
 
            lambda = Temp.Multiply(gi);
  
            // Calcul de P
            Matrix P = new DenseMatrix(N, N);
            Matrix Id = new DenseMatrix(N, N);
 
            for (i = 0; i < N; i++)
            {
                for (j = 0; j < N; j++)
                {
                    Id[i, j] = 0;

                    if (i == j)
                    {
                        Id[i, j] = 1;
                    }
                }
            }

            P = Phi.Inverse().Multiply(Phi2);

            P.Multiply(-0.5 * Math.Pow(sigma, 2));

            Temp = Id.Clone();
            Temp.Multiply(r);

            P.Add(Temp);

            Temp = Phi.Inverse().Multiply(Phi1);
 
            Temp.Multiply(r - 0.5 * Math.Pow(sigma, 2));

            P.Subtract(Temp);

            // Calcul de H et de G
            Matrix H = new DenseMatrix(N, N);
            Matrix G = new DenseMatrix(N, N);

            H = Id.Clone();

            Temp = P.Clone();

            Temp.Multiply(deltaT * (1-theta));

            H.Add(Temp);

            G = Id.Clone();

            Temp = P.Clone();
            Temp.Multiply(deltaT * theta);

            G.Subtract(Temp);
 
            // Boucle de k=1 à k=M
            for (i = 1; i <= M; i++)
            {
                H.Multiply(lambda,lambda);
                Temp = G.Inverse();

                lambda = Temp.Multiply(lambda);
  
                t += deltaT;

                lambda[0] = alpha(t);
                lambda[N - 1] = beta(t, S, r, T, E);
 
                lambda = Phi.Inverse().Multiply(lambda);
            }

            for (i = 0; i < N; i++)
            {
                Prix += lambda[i] * phi(x-ksi[i], c);
            }

            return Prix;

        }

Are there any errors?
Thanks in advance.

 

[Tsotne]

The code seems to be ok at the first glance. Make sure you're supplying arguments not causing the code to result in logic error.