% This is the Matlab code used to create Figures 3, 5, 7, 8 and 9 in "What Are The Odds? Option-Based % Forecasts of FOMC Target Changes" in the November/December issue of the Review, Volume 88(6). % % The program uses MLE to estimate the probabilities that the Fed will choose a particular funds rate target from % a range of targets, using options on Fed Funds Futures contracts and the restrictions that the probabilities % are all greater than or equal to 0, that the sum of the probabilities is one, and that the Futures price % assumes a 1 bp/month risk premium. % % This program requires the defined procedure defined below: % function ssres = ssres(phat) % global X aop; % ssres = (aop - X*phat)'*(aop - X*phat); % % This function must accompany the Matlab code for the estimation procedure to work. % % Input Data are described in the accompanying data file. clear all; global X aop; %--------------------------------------------------- % DEFINE INPUTS, OUTPUTS, AND GENERAL SPECIFICATIONS %--------------------------------------------------- % Input File filename = 'G:\Temp\Matlab\December2006(10.4).xls'; % Output File outputfilename = 'G:\Temp\Matlab\December2006(MatlabOutput)2.xls'; %These are the worksheet names from the input Excel file template worksheetname1 = 'Futures and TBill'; worksheetname2 = 'Bloomberg Notation'; %Replace numerical columns with the Excel column identification characters colindex = char('a','b','c','d','e','f','g','h','i','j','k','l','m','n','o',... 'p','q','r','s','t','u','v','w','x','y','z','aa','ab','ac','ad','ae','af',... 'ag','ah','ai','aj','ak','al','am','an','ao','ap','aq','ar','as','at','au',... 'av','aw','ax','ay','az','ba','bb','bc','bd','be','bf','bg','bh','bi','bj',... 'bk','bl','bm','bn','bo','bp','bq','br','bs','bt','bu','bv','bw','bx','by','bz'); %------------------------------------ % DEFINE CONSTANTS FROM INPUT DATA %------------------------------------ nop = xlsread(char(filename),worksheetname1,'b9'); %This is the total number of options in the file for each date [N,firstday] = xlsread(char(filename),worksheetname1,'a14'); %First date of input dataset %FOMC Dates% lastfomc = xlsread(char(filename),worksheetname1,'b2:d2'); lastfomcmonth = lastfomc(1,1); lastfomcday = lastfomc(1,2); lastfomcyear = lastfomc(1,3); %Contract Expiration Dates% contractmonth = xlsread(char(filename),worksheetname1,'c8'); contractyear = xlsread(char(filename),worksheetname1,'d8'); contractday = eomday(contractyear,contractmonth); %Gives the last day of the month of the current contract %Number of Scenarios and Potential high-end Fed Funds rate targets (We will generally use this program for only one scenario)% numscen = xlsread(char(filename),worksheetname2,'b15'); %Number of Scenarios ffflowhigh = xlsread(char(filename),worksheetname1,strcat('a4:',colindex(numscen,:),'5')); %This gives the highest and lowest considered Fed Funds targets for the forecasted meeting %Define the Dimensions of the Data Set% firstdayserial = datenum(firstday); %Replaces the first date with a numerical counting vector starting at one numdays = xlsread(char(filename),worksheetname1,'D9'); alldata = xlsread(char(filename),worksheetname1,strcat('b13:',colindex(4+nop,:),num2str(13+numdays))); %This is the total # of calls and puts traded on the contract %Number of Options nop_alt = sum(~isnan(alldata(:,4:end)),2); %This gives the total number of options traded for each day. nop_alt1 = ~isnan(alldata(:,4:end)); % 1 if options traded, zero if not. %----------------------------------- % DEFINE OPTIONS AT EACH STRIKE PRICE %----------------------------------- aox = xlsread(char(filename),worksheetname1,strcat('e8:',colindex(4+nop,:),'8'))'; %This reads in the list of potential strike prices. They are in 6.25 basis point intervals. [N,aodc] = xlsread(char(filename),worksheetname1,strcat('e7:',colindex(4+nop,:),'7')); %This is just the string "C" or "P" for the strikes depending on if they are calls or puts aodc = aodc'; ifmissing = sum(isnan(alldata),2)>(nop-5) | isnan(alldata(:,1)); %This calculates the number of missing days in the data [N,alldates] = xlsread(char(filename),worksheetname1,strcat('A','13',':','A',num2str(13+numdays))); %This reads in the number of days from the start date to the end dayexp = datenum(contractyear,contractmonth,contractday) - datenum(alldates); %Number of days left on contract until expiration (ALL DAYS, NOT BUSINESS DAYS!!!) %------------------------------------------------------------------- % NORMALIZE MEASUREMENT DATA TO BE DIRECTLY COMPARABLE FOR ANALYSIS % ------------------------------------------------------------------ %Data modification% fffutd = 100 - alldata(:,1); % 100-Fed Funds Future Daily Closing Price tbillbloom = floor(100*alldata(:,2) + .005)/100; %Round 6-month T-bill rate to the nearest integer (lower-bound). pricebloom = 100 - ((100*dayexp/360).*(tbillbloom/100)); %This calculation takes account of the risk premium tbill = floor(1000000*((100./pricebloom).^(360./dayexp)-1+0.0000005))/1000000; %This is the price of the 6-month T-bill taking account of the interest rate differential d = exp(-(tbill).*(dayexp/360)); datadays = numdays - sum(ifmissing); %Total number of days in sample - the number of missing values %Create zeros vectors for future output% datadatesserial = zeros(datadays,1); fffutdnomiss = zeros(datadays,1); fffutdTPnomiss = zeros(datadays,1); %Loop through the number of scenarios (probably only one)% for iii = 1:numscen tff = (ffflowhigh(1,iii):.25:ffflowhigh(2,iii))'; %Creates an index of potential Fed Fund target rates in 25 basis point increments nfff = size(tff,1); %Dimensions of the potential fed funds rates (one column and variable rows) %If we use the contract that expires in the same month as the FOMC meeting, %then we use the futures-implied Fed Funds rate, weighted by the %number of days in the month prior to the meeting, and the matrix %of potential fed funds rate targets weighted by the number of days %in the month after the target change in order to compute the %options payoff matrix. if contractmonth == lastfomcmonth fffbegmonth = xlsread(char(filename),worksheetname1,'a7'); fff = fffbegmonth*((lastfomcday-1)/contractday) + (1-((lastfomcday-1)/contractday))*tff; else fff = tff; end %---------------------------------------------------------- %DEFINE AND IMPOSE RESTRICTIONS FOR MLE ESTIMATION PROCEDURE %----------------------------------------------------------- %Create empty matrices for output% phatresSum1 = zeros(datadays,nfff); phatresEqfff = zeros(datadays,nfff); phatresTP = zeros(datadays,nfff); phatresSum1_ge0 = zeros(datadays,nfff); phatresEqfff_ge0 = zeros(datadays,nfff); phatresTP_ge0 = zeros(datadays,nfff); daydata = 1; % Counter %Loop through each day in the sample% for jjj = 1:numdays if ifmissing(jjj) == 0 datadatesserial(daydata,1) = datenum(alldates(jjj,1)); % Put today's date in the zeros vector dateswdata(daydata,1) = alldates(jjj,1); op1 = alldata(jjj,4:size(alldata,2))'; % Today's option traded allopdata = alldata(:,4:size(alldata,2)); % All option data op = allopdata(jjj,~isnan(allopdata(jjj,:)))'; aop = op./d(jjj,1); %Option divided by the interest rate differential X = zeros(nop_alt(jjj),size(fff,1)); %X-Matrix: Columns are Fed Funds Rates. Rows are Strike Prices counter = 1; %COMPUTE PAYOFF MATRIX: This section can be confusing. Essentially what we are interested in is %whether or not the price of a call option is less than the futures price. If it is, then we want to exercise %the option. The prices are not computed as 100-price in this data set though. So for a %call option, we want to exercise it only if the strike price is greater than %the futures price, because then 100-price will be lower. The reverse is true for a put. for kkk = 1:nop %Loop through each option at each Futures price if nop_alt1(jjj,kkk) == 1 if strncmp(aodc(kkk,1),'C',1) | strncmp(aodc(kkk,1),'c',1) %For call options X(counter,:) = max(aox(kkk,1)-fff',0); counter = counter + 1; elseif strncmp(aodc(kkk,1),'P',1) | strncmp(aodc(kkk,1),'p',1) %For put options X(counter,:) = max(fff'-aox(kkk,1),0); counter = counter + 1; else error('error with options'); end end end fffutdnomiss(daydata,1) = fffutd(jjj,1); %100-Fed Funds rate fffutdTPnomiss(daydata,1) = fffutd(jjj,1) - (dayexp(jjj)/30)/100; %100-FFF, using risk premium %------------------------------------- % ESTIMATE THE PROBABILITIES USING MLE %------------------------------------- %Restriction 1: Default (probs sum to 1)% r1 = zeros(size(X,2),1); r1(size(X,2),1) = 1; r2 = zeros(size(X,2),size(X,2)-1); r2(1:size(X,2)-1,:) = eye(size(X,2)-1); r2(size(X,2),:) = -1*ones(1,size(X,2)-1); phatresSum1(daydata,:) = (r1 + r2*(inv((X*r2)'*(X*r2))*(X*r2)'*(aop-X*r1)))'; %Estimates %Restriction 2: Each individual probability is greater than or equal to 0% Aeq = ones(1,size(X,2)); beq = 1; A = -eye(size(X,2)); b = zeros(size(X,2),1); phatresSum1_ge0(daydata,:) = fmincon(@ssres,phatresSum1(daydata,:)',A,b,Aeq,beq,[],[],[],optimset('TolFun',1e-12,'TolCon',1e-12,'Display','off'))'; %Estimates %Restriction 3: mean PDF = futures price% D = eye(size(fff,1)); D(size(fff,1)-1,:)=ones(1,size(fff,1)); D(size(fff,1),:)=fff'; d1 = zeros(size(fff,1),1); d1(size(fff,1)-1,1)=1; d1(size(fff,1),1)=fffutd(jjj,1); d2 = zeros(size(fff,1),size(fff,1)-2); d2(1:size(fff,1)-2,1:size(fff,1)-2) = eye(size(fff,1)-2); r1 = inv(D)*d1; r2 = inv(D)*d2; phatresEqfff(daydata,:) = (r1 + r2*(inv((X*r2)'*(X*r2))*(X*r2)'*(aop-X*r1)))'; %Estimates %Restriction 4: Probability are grater than or equal to 0 AND mean PDF = futures price Aeq = ones(2,size(X,2)); Aeq(2,:) = fff'; beq = [1 fffutd(jjj,1)]'; phatresEqfff_ge0(daydata,:) = fmincon(@ssres,phatresEqfff(daydata,:)',A,b,Aeq,beq,[],[],[],optimset('TolFun',1e-12,'TolCon',1e-12,'Display','off'))'; %Estimates %Restriction 5: Mean PDF = premium adjusted futures price% d1(size(fff,1),1)=fffutd(jjj,1) - (dayexp(jjj)/30)/100; r1 = inv(D)*d1; phatresTP(daydata,:) = (r1 + r2*(inv((X*r2)'*(X*r2))*(X*r2)'*(aop-X*r1)))'; %Estimates %Restriction 6: Probabilities are greater than or equal to 0 AND mean PDF = premium adjusted futures price% beq = [1 fffutd(jjj,1) - (dayexp(jjj)/30)/100]'; %Adjust for risk premium phatresTP_ge0(daydata,:) = fmincon(@ssres,phatresTP(daydata,:)',A,b,Aeq,beq,[],[],[],optimset('TolFun',1e-12,'TolCon',1e-12,'Display','off'))'; %Estimates daydata = daydata+1; end %This ends the midding date if statement end %This ends the jjj loop %------------------------------- % OUTPUT AND GRAPH RESULTS %------------------------------- titnames = char('default (probs sum to 1)','mean PDF = futures price','mean PDF = premium adjusted futures','probs ge 0','probs ge 0, mean PDF = futures price','probs ge 0, mean PDF = premium adjusted futures'); %Empty matrices for each of the restriction results% phats = zeros(6,size(phatresSum1,1),size(phatresSum1,2)); phats(1,:,:)=phatresSum1; phats(2,:,:)=phatresEqfff; phats(3,:,:)=phatresTP; phats(4,:,:)=phatresSum1_ge0; phats(5,:,:)=phatresEqfff_ge0; phats(6,:,:)=phatresTP_ge0; %Loop through each scenario to print output and graph results% figure, orient landscape; for jjj = 1:6 %Different loop for each restriction subplot(2,3,jjj); plot(datadatesserial,squeeze(phats(jjj,:,:))); xlim([datadatesserial(1,1) datadatesserial(size(datadatesserial,1),1)]); ylim([-.2 1.2]); set(gca,'XTick',datadatesserial(1,1):21:datadatesserial(size(datadatesserial,1),1)+5); set(gca,'YTick',-.2:.1:1.2); datetick('x',6,'keepticks','keeplimits'); legend(num2str((ffflowhigh(1,iii):.25:ffflowhigh(2,iii))'),'Location','Best'); legend('boxoff'); title(titnames(jjj,:)); filenamep = char(outputfilename); worksheetnamep = strcat('scenario ',num2str(jjj)); temptff = tff'/100; expprobs = squeeze(phats(jjj,:,:)); expectation = expprobs*temptff'; expected = cellstr('Expectation'); xlswrite(filenamep,temptff,worksheetnamep,'b1'); %Output potential target Fed Funds Rates xlswrite(filenamep,dateswdata,worksheetnamep,'a2'); %Output Date Vector (Column 1) xlswrite(filenamep,squeeze(phats(jjj,:,:)),worksheetnamep,'b2'); xlswrite(filenamep,expectation,worksheetnamep,strcat((colindex(size(temptff,2)+2)),'2')); xlswrite(filenamep,expected,worksheetnamep,strcat((colindex(size(temptff,2)+2)),'1')); end %This ends the jjj loop end %This ends the scenario loop