// d B(r0,lambda0,t)/dt = mean_reversion_r(r0-r)*dB/dr + 0.5*sigma_r^2*d^2B/dr^2 + // 0.5*sigma_lambda^2*d^2B/dlambda^2 + rho*sigma_r*sigma_lambda*d^2B/dr*dlambda // -(r+lambda)*B // g(t,r,lambda)*lambda + f(t,r,lambda) // where g(t,r,lambda) - expected value of integral // f(t,r,lambda) - continuouse source/sink // // //This code assumes relatively large diffusion coefficient (I used 0.3) //if you make it small, than well-known oscillatory effect arise //due to central difference approximartion of the first //order derivative, so keep this in mind and if you want small vall, than //change approximation of the first order derivative //or use "particle methods" //also you need "dumka.c" to run this //depending on your compiler you can attach it in different manner //the code suppose tobeself-explanatory,rather than very efficient, but if //you can sqeese more than 20% efficiency out of it-let me know please. #include "stdafx.h"//commen this line out if you are not unver Miceosoft Visual Studio #include #include #ifndef max #define max(a,b) #endif FILE *payoffFunction; FILE *solution; FILE *solution_exact; double t_old=-0.1e0; int nPointsX=200; int nPointsY=200; double x_min=-0.0; double x_max=1.0; double l_min=-0.0; double l_max=1.0; double intervalX=x_max-x_min; double intervalL=l_max-l_min; double dx=intervalX/(double) (nPointsX-1); double dl=intervalL/(double) (nPointsY-1); long counter_call_back=0; double mean_Reverting_r=0.05e0; double mean_Reverting_l=0.05e0; double rho = -0.2; double r0=0.09e0;//0.09e0; double rate_of_meanreversion_r=2.0; double rate_of_meanreversion_l=2.0; double *x_arr, *l_arr; double pi=acos(-1.0e0); double final_payoff(double T, double r, double l) { return 1.0e0; } double continuouse_payoff(double t, double r, double l) { return 0.0e0; } double default_payoff(double t, double r, double l, double pi) { return 0.0e0; } double expected_default_payoff(double t, double r, double l) { return 0.0e0; } double sigma_r(double t, double x){ double sigma=0.3e0; return sigma; } double sigma_l(double t, double l){ double sigma=0.3e0; return sigma; } void rhs(long neqn, double* y, double t,double *result){ //double y_Lft, y_Rgt, y_Top, y_Bottom, y_Center; double sigma_r_tmp, sigma_l_tmp; double drift_r=0.0e0; double drift_l=0.0e0; double dx_sq_inv=1.0e0/(dx*dx); double dl_sq_inv=1.0e0/(dl*dl); double dx_inv=1.0e0/dx; double dl_inv=1.0e0/dl; double dx_dl_inv=dx_inv*dl_inv; //--------------- not finished boundary conditions yet ----------------- int i,j,k; i=0; for ( j=0; j