| //HJM.cpp |
| //Routine to setup HJM framework. |
| //Authors: Mark Broadie, Jatin Dewanwala, Columbia University |
| //Collaborator: Mikhail Smelyanskiy, Intel |
| //Based on hjm_simn.xls created by Mark Broadie |
| |
| #include <stdio.h> |
| #include <stdlib.h> |
| #include <math.h> |
| #include <time.h> |
| |
| #include "nr_routines.h" |
| #include "HJM.h" |
| #include "HJM_type.h" |
| |
| |
| |
| int HJM_SimPath_Yield(FTYPE **ppdHJMPath, int iN, int iFactors, FTYPE dYears, FTYPE *pdYield, FTYPE **ppdFactors, |
| long *lRndSeed); |
| int HJM_SimPath_Forward(FTYPE **ppdHJMPath, int iN, int iFactors, FTYPE dYears, FTYPE *pdForward, FTYPE *pdTotalDrift, |
| FTYPE **ppdFactors, long *lRndSeed); |
| int HJM_Yield_to_Forward(FTYPE *pdForward, int iN, FTYPE *pdYield); |
| int HJM_Factors(FTYPE **ppdFactors,int iN, int iFactors, FTYPE *pdVol, FTYPE **ppdFacBreak); |
| int HJM_Drifts(FTYPE *pdTotalDrift, FTYPE **ppdDrifts, int iN, int iFactors, FTYPE dYears, FTYPE **ppdFactors); |
| int HJM_Correlations(FTYPE **ppdHJMCorr, int iN, int iFactors, FTYPE **ppdFactors); |
| int HJM_Forward_to_Yield(FTYPE *pdYield, int iN, FTYPE *pdForward); |
| int Discount_Factors(FTYPE *pdDiscountFactors, int iN, FTYPE dYears, FTYPE *pdRatePath); |
| //int Discount_Factors_early_exit(FTYPE *pdDiscountFactors, int iN, FTYPE dYears, FTYPE *pdRatePath, int iSwapStartTimeIndex); |
| |
| int HJM_SimPath_Yield(FTYPE **ppdHJMPath, //Matrix that stores generated HJM path (Output) |
| int iN, //Number of time-steps |
| int iFactors, //Number of factors in the HJM framework |
| FTYPE dYears, //Number of years |
| FTYPE *pdYield, //Input yield curve (at t=0) for dYears (iN time steps) |
| FTYPE **ppdFactors, //Matrix of Factor Volatilies |
| long *lRndSeed) |
| { |
| //This function returns a single generated HJM Path for the given inputs |
| |
| int iSuccess = 0; //return variable |
| |
| FTYPE *pdForward; //Vector that will store forward curve computed from given yield curve |
| FTYPE **ppdDrifts; //Matrix that will store drifts for different maturities for each factor |
| FTYPE *pdTotalDrift; //Vector that stores total drift for each maturity |
| |
| pdForward = dvector(0, iN-1); |
| ppdDrifts = dmatrix(0, iFactors-1, 0, iN-2); |
| pdTotalDrift = dvector(0, iN-2); |
| |
| //generating forward curve at t=0 from supplied yield curve |
| iSuccess = HJM_Yield_to_Forward(pdForward, iN, pdYield); |
| if (iSuccess!=1) |
| { |
| free_dvector(pdForward, 0, iN-1); |
| free_dmatrix(ppdDrifts, 0, iFactors-1, 0, iN-2); |
| free_dvector(pdTotalDrift, 0, iN-1); |
| return iSuccess; |
| } |
| |
| //computation of drifts from factor volatilities |
| iSuccess = HJM_Drifts(pdTotalDrift, ppdDrifts, iN, iFactors, dYears, ppdFactors); |
| if (iSuccess!=1) |
| { |
| free_dvector(pdForward, 0, iN-1); |
| free_dmatrix(ppdDrifts, 0, iFactors-1, 0, iN-2); |
| free_dvector(pdTotalDrift, 0, iN-1); |
| return iSuccess; |
| } |
| |
| //generating HJM Path |
| iSuccess = HJM_SimPath_Forward(ppdHJMPath, iN, iFactors, dYears, pdForward, pdTotalDrift,ppdFactors, lRndSeed); |
| if (iSuccess!=1) |
| { |
| free_dvector(pdForward, 0, iN-1); |
| free_dmatrix(ppdDrifts, 0, iFactors-1, 0, iN-2); |
| free_dvector(pdTotalDrift, 0, iN-1); |
| return iSuccess; |
| } |
| |
| free_dvector(pdForward, 0, iN-1); |
| free_dmatrix(ppdDrifts, 0, iFactors-1, 0, iN-2); |
| free_dvector(pdTotalDrift, 0, iN-1); |
| iSuccess = 1; |
| return iSuccess; |
| } |
| |
| |
| int HJM_Yield_to_Forward (FTYPE *pdForward, //Forward curve to be outputted |
| int iN, //Number of time-steps |
| FTYPE *pdYield) //Input yield curve |
| { |
| //This function computes forward rates from supplied yield rates. |
| |
| int iSuccess=0; |
| int i; |
| |
| //forward curve computation |
| pdForward[0] = pdYield[0]; |
| for(i=1;i<=iN-1; ++i){ |
| pdForward[i] = (i+1)*pdYield[i] - i*pdYield[i-1]; //as per formula |
| //printf("pdForward: %f = (%d+1)*%f - %d*%f \n", pdForward[i], i, pdYield[i], i, pdYield[i-1]); |
| } |
| iSuccess=1; |
| return iSuccess; |
| } |
| |
| |
| int HJM_Factors(FTYPE **ppdFactors, //Output matrix that stores factor volatilities for different maturities |
| int iN, |
| int iFactors, |
| FTYPE *pdVol, //Input vector of total volatilities for different maturities |
| FTYPE **ppdFacBreak) //Input matrix of factor weights for each maturity |
| { |
| //This function computes individual volatilities for each factor for different maturities. |
| //The function is called when the user inputs total volatility data and the weight distribution |
| //according to which the total variance has to be split accross various factors. |
| |
| //For instance, the user may supply Maturity: 1 2 3 4 |
| //total vol (pdVol) as Sigma: 1.35%, 1.30%, 1.25%, 1.20%,.... |
| //and the weight breakdown (ppdFacBreak) as Factor 1: 0.55, 0.60, 0.65, 0.69,.... |
| // Factor 2: 0.44, 0.39, 0.34, 0.30,.... |
| // Factor 3: 0.01, 0.01, 0.01, 0.01,.... |
| //Note that the weights add up to 1 in each case. Also, the weights are based on variance not volatility. |
| |
| //Based on these inputs, the function will calculate individual volatilties for each factor for each maturity. |
| //The output (ppdFactors) may look something like: Maturity: 1 2 3 4 |
| // Factor 1: 1.00% 1.00% 1.00% 1.00% |
| // Factor 2: 0.90% 0.82% 0.74% 0.67% |
| // Factor 3: 0.10% 0.08% 0.05% 0.03% |
| // (Please note that in this example the figures have been rounded and therefore may not be exact.) |
| |
| int i,j; //looping variables |
| int iSuccess = 0; |
| |
| //Computation of factor volatilities |
| for(i = 0; i<=iFactors-1; ++i) |
| for(j=0; j<=iN-2;++j) |
| ppdFactors[i][j] = sqrt((ppdFacBreak[i][j])*(pdVol[j])*(pdVol[j])); |
| |
| iSuccess =1; |
| return iSuccess; |
| } |
| |
| |
| int HJM_Drifts(FTYPE *pdTotalDrift, //Output vector that stores the total drift correction for each maturity |
| FTYPE **ppdDrifts, //Output matrix that stores drift correction for each factor for each maturity |
| int iN, |
| int iFactors, |
| FTYPE dYears, |
| FTYPE **ppdFactors) //Input factor volatilities |
| { |
| //This function computes drift corrections required for each factor for each maturity based on given factor volatilities |
| |
| int iSuccess =0; |
| int i, j, l; //looping variables |
| FTYPE ddelt = (FTYPE) (dYears/iN); |
| FTYPE dSumVol; |
| |
| //computation of factor drifts for shortest maturity |
| for (i=0; i<=iFactors-1; ++i) |
| ppdDrifts[i][0] = 0.5*ddelt*(ppdFactors[i][0])*(ppdFactors[i][0]); |
| |
| //computation of factor drifts for other maturities |
| for (i=0; i<=iFactors-1;++i) |
| for (j=1; j<=iN-2; ++j) |
| { |
| ppdDrifts[i][j] = 0; |
| for(l=0;l<=j-1;++l) |
| ppdDrifts[i][j] -= ppdDrifts[i][l]; |
| dSumVol=0; |
| for(l=0;l<=j;++l) |
| dSumVol += ppdFactors[i][l]; |
| ppdDrifts[i][j] += 0.5*ddelt*(dSumVol)*(dSumVol); |
| } |
| |
| //computation of total drifts for all maturities |
| for(i=0;i<=iN-2;++i) |
| { |
| pdTotalDrift[i]=0; |
| for(j=0;j<=iFactors-1;++j) |
| pdTotalDrift[i]+= ppdDrifts[j][i]; |
| } |
| |
| iSuccess=1; |
| return iSuccess; |
| } |
| |
| int HJM_SimPath_Forward(FTYPE **ppdHJMPath, //Matrix that stores generated HJM path (Output) |
| int iN, //Number of time-steps |
| int iFactors, //Number of factors in the HJM framework |
| FTYPE dYears, //Number of years |
| FTYPE *pdForward, //t=0 Forward curve |
| FTYPE *pdTotalDrift, //Vector containing total drift corrections for different maturities |
| FTYPE **ppdFactors, //Factor volatilities |
| long *lRndSeed) //Random number seed |
| { |
| //This function computes and stores an HJM Path for given inputs |
| |
| int iSuccess = 0; |
| int i,j,l; //looping variables |
| |
| FTYPE ddelt; //length of time steps |
| FTYPE dTotalShock; //total shock by which the forward curve is hit at (t, T-t) |
| FTYPE *pdZ; //vector to store random normals |
| |
| ddelt = (FTYPE)(dYears/iN); |
| |
| pdZ = dvector(0, iFactors -1); //assigning memory |
| |
| for(i=0;i<=iN-1;++i) |
| for(j=0;j<=iN-1;++j) |
| ppdHJMPath[i][j]=0; //initializing HJMPath to zero |
| |
| //t=0 forward curve stored iN first row of ppdHJMPath |
| for(i=0;i<=iN-1; ++i) |
| ppdHJMPath[0][i] = pdForward[i]; |
| |
| //Generation of HJM Path |
| for (j=1;j<=iN-1;++j) |
| { |
| |
| for (l=0;l<=iFactors-1;++l) |
| pdZ[l]= CumNormalInv(RanUnif(lRndSeed)); //shocks to hit various factors for forward curve at t |
| |
| for (l=0;l<=iN-(j+1);++l) |
| { |
| dTotalShock = 0; |
| for (i=0;i<=iFactors-1;++i) |
| dTotalShock += ppdFactors[i][l]* pdZ[i]; |
| ppdHJMPath[j][l] = ppdHJMPath[j-1][l+1]+ pdTotalDrift[l]*ddelt + sqrt(ddelt)*dTotalShock; |
| //as per formula |
| } |
| } |
| |
| free_dvector(pdZ, 0, iFactors -1); |
| iSuccess = 1; |
| return iSuccess; |
| } |
| |
| |
| |
| |
| int HJM_Correlations(FTYPE **ppdHJMCorr,//Matrix that stores correlations among factor volatilities for different maturities |
| int iN, |
| int iFactors, |
| FTYPE **ppdFactors) |
| { |
| //This function is based on factor.xls created by Mark Broadie |
| //The function computes correlations between factor volatilities for different maturities |
| |
| int iSuccess = 0; |
| int i, j, l; //looping variables |
| FTYPE *pdTotalVol; //vector that stores total volatility data for different maturities |
| FTYPE **ppdWeights; //matrix that stores ratio of each factor to total volatility for different maturities |
| |
| pdTotalVol = dvector(0,iN-2); |
| ppdWeights = dmatrix(0, iFactors-1,0, iN-2); |
| |
| //Total Volatility computed from given factor volatilities |
| for(i=0;i<=iN-2;++i) |
| { |
| pdTotalVol[i]=0; |
| for(j=0;j<=iFactors-1;++j) |
| pdTotalVol[i] += ppdFactors[j][i]*ppdFactors[j][i]; |
| pdTotalVol[i] = sqrt(pdTotalVol[i]); |
| } |
| |
| //Weights computed |
| for(i=0;i<=iN-2;++i) |
| for(j=0;j<=iFactors-1;++j) |
| ppdWeights[j][i] = ppdFactors[j][i]/pdTotalVol[i]; |
| |
| //Output matrix initialized to zero |
| for(i=0;i<=iN-2;++i) |
| for(j=0;j<=iN-2;++j) |
| ppdHJMCorr[i][j]=0; |
| |
| //Correlations computed |
| for(i=0;i<=iN-2;++i) |
| for(j=i;j<=iN-2;++j) |
| for(l=0;l<=iFactors-1;++l) |
| ppdHJMCorr[i][j] += ppdWeights[l][i]*ppdWeights[l][j]; |
| |
| free_dvector(pdTotalVol, 0,iN-2); |
| free_dmatrix(ppdWeights, 0, iFactors-1,0, iN-2); |
| iSuccess = 1; |
| return iSuccess; |
| } |
| |
| |
| int HJM_Forward_to_Yield (FTYPE *pdYield, //Output yield curve |
| int iN, |
| FTYPE *pdForward) //Input forward curve |
| { |
| //This function computes yield rates from supplied forward rates. |
| |
| int iSuccess=0; |
| int i; |
| |
| //t=0 yield curve |
| pdYield[0] = pdForward[0]; |
| for(i=1;i<=iN-1; ++i) |
| pdYield[i] = (i*pdYield[i-1] + pdForward[i])/(i+1); |
| |
| iSuccess=1; |
| return iSuccess; |
| } |
| |
| int Discount_Factors(FTYPE *pdDiscountFactors, |
| int iN, |
| FTYPE dYears, |
| FTYPE *pdRatePath) |
| { |
| int i,j; //looping variables |
| int iSuccess; //return variable |
| |
| FTYPE ddelt; //HJM time-step length |
| ddelt = (FTYPE) (dYears/iN); |
| |
| //initializing the discount factor vector |
| for (i=0; i<=iN-1; ++i) |
| pdDiscountFactors[i] = 1.0; |
| |
| for (i=1; i<=iN-1; ++i) |
| for (j=0; j<=i-1; ++j) |
| pdDiscountFactors[i] *= exp(-pdRatePath[j]*ddelt); |
| |
| iSuccess = 1; |
| return iSuccess; |
| } |
| |
| int Discount_Factors_opt(FTYPE *pdDiscountFactors, |
| int iN, |
| FTYPE dYears, |
| FTYPE *pdRatePath) |
| { |
| int i,j; //looping variables |
| int iSuccess; //return variable |
| |
| FTYPE ddelt; //HJM time-step length |
| ddelt = (FTYPE) (dYears/iN); |
| |
| FTYPE *pdexpRes; |
| pdexpRes = dvector(0,iN-2); |
| |
| //initializing the discount factor vector |
| for (i=0; i<=iN-1; ++i) |
| pdDiscountFactors[i] = 1.0; |
| |
| //precompute the exponientials |
| for (j=0; j<=(i-2); ++j){ pdexpRes[j] = -pdRatePath[j]*ddelt; } |
| for (j=0; j<=(i-2); ++j){ pdexpRes[j] = exp(pdexpRes[j]); } |
| |
| for (i=1; i<=iN-1; ++i) |
| for (j=0; j<=i-1; ++j) |
| pdDiscountFactors[i] *= pdexpRes[j]; |
| |
| free_dvector(pdexpRes, 0, iN-2); |
| iSuccess = 1; |
| return iSuccess; |
| } |
| |
| |
| // *********************************************************************** |
| // *********************************************************************** |
| // *********************************************************************** |
| int Discount_Factors_Blocking(FTYPE *pdDiscountFactors, |
| int iN, |
| FTYPE dYears, |
| FTYPE *pdRatePath, |
| int BLOCKSIZE) |
| { |
| int i,j,b; //looping variables |
| int iSuccess; //return variable |
| |
| FTYPE ddelt; //HJM time-step length |
| ddelt = (FTYPE) (dYears/iN); |
| |
| FTYPE *pdexpRes; |
| pdexpRes = dvector(0,(iN-1)*BLOCKSIZE-1); |
| //precompute the exponientials |
| for (j=0; j<=(iN-1)*BLOCKSIZE-1; ++j){ pdexpRes[j] = -pdRatePath[j]*ddelt; } |
| for (j=0; j<=(iN-1)*BLOCKSIZE-1; ++j){ pdexpRes[j] = exp(pdexpRes[j]); } |
| |
| |
| //initializing the discount factor vector |
| for (i=0; i<(iN)*BLOCKSIZE; ++i) |
| pdDiscountFactors[i] = 1.0; |
| |
| for (i=1; i<=iN-1; ++i){ |
| //printf("\nVisiting timestep %d : ",i); |
| for (b=0; b<BLOCKSIZE; b++){ |
| //printf("\n"); |
| for (j=0; j<=i-1; ++j){ |
| pdDiscountFactors[i*BLOCKSIZE + b] *= pdexpRes[j*BLOCKSIZE + b]; |
| //printf("(%f) ",pdexpRes[j*BLOCKSIZE + b]); |
| } |
| } // end Block loop |
| } |
| |
| free_dvector(pdexpRes, 0,(iN-1)*BLOCKSIZE-1); |
| iSuccess = 1; |
| return iSuccess; |
| } |
| |
| |
| |