Message-ID: <2280763.1075863438101.JavaMail.evans@thyme>
Date: Thu, 16 Aug 2001 16:16:37 -0700 (PDT)
From: j.kaminski@enron.com
To: vkaminski@aol.com
Subject: FW: HJM forward curve modeling
Mime-Version: 1.0
Content-Type: text/plain; charset=ANSI_X3.4-1968
Content-Transfer-Encoding: 7bit
X-From: Kaminski, Vince J </O=ENRON/OU=NA/CN=RECIPIENTS/CN=VKAMINS>
X-To: 'vkaminski@aol.com'
X-cc: 
X-bcc: 
X-Folder: \VKAMINS (Non-Privileged)\Kaminski, Vince J\Sent Items
X-Origin: Kaminski-V
X-FileName: VKAMINS (Non-Privileged).pst



 -----Original Message-----
From: 	Lu, Zimin  
Sent:	Thursday, August 16, 2001 9:58 AM
To:	Halliburton, Tom
Subject:	HJM forward curve modeling




Tom, 

The attached C routine is a HJM forward curve modeling module in the storage model.  
Please study it and then we get together to talk about how to use it for your optimization
model. Pay attention to the blue section. 


Zimin


-------------------

void HJMsim( int nRun, int nfactor, int toEndMonth, \
			 double *F, double *vol, double s0, double *svol, double *correl, 
			 double *LHP, double *mean_spread,
			 int freq,int freq1, int *bin_dim, double **bin, int bin_dim2, \
 			 double **a, double **c, double **d, float ***prob, int *index0){

long traj;
int  i, j, k, it,is, is2,index,month, *pBin, Nt, index_old;
double  *price, **fwdCv, **fwdCv2, *ffVol, *bin2, temp1,temp2, DT, sqrtDT;
double correlp, x;
float  ***nrv;
double F0_0,F0, spot0,spot1, spot, prompt,prompt2, spread; // two state variables
double **L,**LH, **Correl;
FILE   *fp;
double var_s, var_F0, delta, b,cc;
int    nperiods;

if(RECORD){
   fp=fopen("C:\\temp\\curve.txt","w");
   fprintf(fp,"traj\tmonth\tperiod\tspot\t\tprompt\t\tprompt+1\tzeroFwd\n");
}



Nt = freq1 + toEndMonth*freq;


DT=1.0/(12.0*freq);  //freq=4 is weekly step, freq=30 is daily step;
sqrtDT=sqrt(DT);



ffVol =dvector(1,toEndMonth);
LH    =dmatrix(1,toEndMonth, 1,nfactor); // historical factor loadings
Correl=dmatrix(1,toEndMonth, 1,toEndMonth);
L     =dmatrix(1,toEndMonth, 1,nfactor); // current factor loadings


for(i=1;i<=toEndMonth;i++)
	for(j=1;j<=nfactor;j++)
	    LH[i][j]= LHP[j+(i-1)*nfactor];

hisCorrel(toEndMonth, nfactor, LH, Correl);
sffVolatility(toEndMonth, vol, ffVol);	
preHJM(toEndMonth,nfactor,ffVol,Correl,L);

free_dmatrix(LH,1,toEndMonth, 1,nfactor);
free_dmatrix(Correl,1,toEndMonth, 1,toEndMonth);



fwdCv=dmatrix(1,toEndMonth,1, freq);
fwdCv2=dmatrix(1,toEndMonth,1, freq);
price=dvector(0,toEndMonth);
nrv =f3tensor(0,toEndMonth, 1, freq, 1, nfactor);
bin2=dvector(1,bin_dim2);
pBin=ivector(0, Nt);	                //Ns X Nt



//Time=0

    spot=s0;               // bining spot
	prompt=F[1];		//bining prompt
	spread=prompt-spot;

	is=hist3(spot, bin_dim[0],  bin[0]);

	if(is<1) is=1;

	for(i=1;i<=bin_dim2;i++)
		if(is-bin_dim2/2+i-1<1)
			bin2[i]=bin[0][1]-bin[0][is];
		else if(is-bin_dim2/2+i-1> bin_dim[0])
		    bin2[i]=bin[0][bin_dim[0]]-bin[0][is];
		else
			bin2[i]=bin[0][is-bin_dim2/2+i-1]-bin[0][is];   //chose bin_dim2 odd

  

	
	is2 = hist3(spread,bin_dim2, bin2);
	if(is2<1) is2=1;

    index= (is-1)*bin_dim2+is2;

	pBin[0]=index;
	*index0=index;



for(traj=1;traj<=nRun; traj++){

//Time=0	   
    spot=s0;
	price[0]=s0;
    for(i=1;i<=toEndMonth;i++)
		price[i]=F[i];


//	a[pBin[0]][0] +=1.0;
//	c[pBin[0]][0] += spot;
//	d[pBin[0]][0] += F[1];

	a[*index0][0] +=1.0;
	c[*index0][0] += spot;
	d[*index0][0] += F[1];

// forward induction
//n-factor model for forward curve

	for(month=0; month<=toEndMonth; month++){
		if(month!=0)
			nperiods=freq;
		else
			nperiods=freq1;

		for(j=1;j<=nperiods;j++)
			for(k=1;k<=nfactor;k++)
		       nrv[month][j][k]=gasdev(&idum);
	}

//propagate to maturity
	for(month=1; month<=toEndMonth; month++)
		for(it=1;it<=month;it++){
			if(it==1)
				nperiods=freq1;
			else
				nperiods=freq;

			for(j=1;j<=nperiods;j++){//do freq times simulation within a month
			      temp1=ffVol[month-it+1]*ffVol[month-it+1]*0.5*DT;
			      for(temp2=0.0,k=1;k<=nfactor;k++)
			      temp2 +=L[month-it+1][k]*sqrtDT*nrv[it][j][k];
			      price[month] *= exp(-temp1+temp2);
			      fwdCv[month][j]=price[month]; //remember only the prompt month contract evolution
                  if(it==month-1 && month !=1)
				    fwdCv2[it][j]=price[month];//prompt+1, NOTE: maturity month=1 can never be prompt+1 
			}
		
		}//it, month


// month =0 current spot

        it=0;  // initialization

		for(month=0;month<toEndMonth;month++){ //spot price and 2-state binning


            correlp=sqrt(1.0-correl[month+1]*correl[month+1]); // correlation between prompt and seasonal spread
                                                               //correl[1] is shifted for correl[0], current spot vs prompt
		    prompt=fwdCv[month+1][1];

			if(month==toEndMonth-1) //
			    prompt2=prompt;  // end point
			else
                prompt2=fwdCv2[month+1][1];

			//F0_0=2.*prompt-prompt2;
            F0_0=price[month];  // convergency theory

			spot0 = F0_0 + mean_spread[month+1];  // mean_spread  need to declared !!!

            if(month!=0){
				var_s=spot0*spot0*(exp(svol[month+1]/12.0)-1.);
				var_F0=F0_0*F0_0*(exp(ffVol[1]/12.0)-1.);
			}
			else{
				var_s=spot0*spot0*(exp(svol[month+1]*freq1*DT)-1.);
				var_F0=F0_0*F0_0*(exp(ffVol[1]*freq1*DT)-1.);
			}


			b=2.*F0_0*correl[month+1]*ffVol[1];
			if(month!=0)
			  cc=(var_F0-var_s)*12.0;
			else
			  cc=(var_F0-var_s)/(freq1*DT);

			temp1=b*b-4.0*cc;
			if(temp1>TINY)
			  delta = (-b+sqrt(temp1))/2.0;
			else
			  delta = TINY;

			spot1 = mean_spread[month+1];// shifted here due to input range

			if(month!=0)
				nperiods=freq;
			else
				nperiods=freq1;


			for(j=1;j<=nperiods;j++){

				if((delta-1.1*TINY)>0.0){
				   x=correl[month+1]*nrv[month][j][1]+correlp*gasdev(&idum); //correl[month] should start from zero, but shifted here
				   temp2 = delta*sqrtDT*x; 
				   spot1 += temp2;
				}


				prompt=fwdCv[month+1][j];
				if(month==toEndMonth-1)  //spot month
			    	prompt2=prompt;  // end point
				else
                    prompt2=fwdCv2[month+1][j];
               
			
				
				F0= prompt - (nperiods-j)*(prompt2-prompt)/(double)freq;
			
				spot = F0 + spot1; // new price process 11/99

				spread=prompt-spot;

                if(RECORD && month<toEndMonth && traj<=10) //Since toEndMonth itself is not in the deal term because we did toEndMonth++ in main routine
					fprintf(fp,"%d\t%d\t%d\t%f\t%f\t%f\t%f \n", traj,month,j, spot,prompt, prompt2,F0);

			
				it++;
				is=hist3(spot, bin_dim[it],  bin[it]);  ///
				if(is<1) is=1;

				for(i=1;i<=bin_dim2;i++)
					if(is-bin_dim2/2+i-1<1)
						bin2[i]=bin[it][1]-bin[it][is];
					else if(is-bin_dim2/2+i-1> bin_dim[it])
					    bin2[i]=bin[it][bin_dim[it]]-bin[it][is];
					else
						bin2[i]=bin[it][is-bin_dim2/2+i-1]-bin[it][is];   //chose bin_dim2 odd

				
				is2   =hist3(spread,bin_dim2, bin2);
				if(is2<1) is2=1;

				index= (is-1)*bin_dim2+is2;

				pBin[it]=index;
				a[index][it] +=1.0;
				c[index][it] += spot;
				d[index][it] += prompt;

/* = populate the transition matrix = */
  
				prob[pBin[it-1]][pBin[it]][it-1] +=1.0;	//[1..Ns]2 [0..Nt-1]?				??			}//end of j?		}//end of month???} // end of traj??/*  =============== free space ===============  */??if(RECORD)?   fclose (fp);??free_ivector(pBin,0,Nt);?free_dmatrix(fwdCv, 1,toEndMonth,1, freq);?free_dmatrix(fwdCv2, 1,toEndMonth,1, freq);?free_dvector(price,0,toEndMonth);?free_f3tensor(nrv,0,toEndMonth, 1, freq, 1, nfactor);?free_dvector(bin2, 1,bin_dim2);?free_dvector(ffVol,1,toEndMonth);??free_dmatrix(L,1,toEndMonth, 1,nfactor); ????}?/*?	for(i=1;i<nCurvePts;i++)?		for(k=Krange[i];k<=Krange[i+1];k++)?			volume[k]=volumeBin[i]+(k-Krange[i])*injectionRate[i];?*/???	return 1;?}?