
% analyze hedge fund index returns
% 

clear all;
close all; 
addpath('..\Data'); 
addpath('..\Programs');
maxdate = 2014.5 ; % axis for graphs

data_source = 1; % 1 = cfsb; 2 = hfr. Using hfr will require some fixing
create_data = 1; % 1 to create factor data from underlying sources. 0 to load data from text file. 
graphs = 0; % 0 can speed up run by not making all the graphs

% ************
% Create data
% This is turned off -- you just load data as below
% skip to "problem set starts" below
% ************

 % new data from CFSB
 
 if data_source == 1;
    x = load('hedgeindexdotcom_data.csv');
    x = x(2:end,:); % hedge data goes through january, factors stop in december
    T = size(x,1);
    x = x(T:-1:1,:); % data is in reverse order 
    HFdates = x(:,1:3); 
    r = x(:,4:end-1);  % ignore multi-strategy rather than program up a separate sample period. 
    THF = size(x,1); 
    %mo,day,year,Hedge Fund Index,Hedge Fund Index Convertible Arbitrage,Hedge Fund Index Dedicated Short Bias,Hedge Fund Index Emerging Markets,Hedge Fund Index Equity Market Neutral,Hedge Fund Index Event Driven,Hedge Fund Index Distressed,Hedge Fund Index Multi-Strategy,Hedge Fund Index Risk Arbitrage,Hedge Fund Index Fixed Income Arbitrage,Hedge Fund Index Global Macro,Hedge Fund Index Long/Short Equity,Hedge Fund Index Managed Futures,Hedge Fund Index Multi-Strategy
    lbls = {'HFIndex','ConvArb','ShortBias','EmergMkt','EquitMktNeut','EventDriven','Distress','Multi-Strat','RiskArb','BondArb','GlobalMacro','LongShtEqty','MgdFuture'}; %Multi-Strategy
    %           1            2           3         4                5                6              7           8         9           10              11          12                13             14          
    plotdate = x(:,1)/12+x(:,3);   
    plotdate_nav = [plotdate(1)-1/12; plotdate]; %??
    subset = [1 2 3 4 5 6 7 8 9 10 11 12];
    allindx = 1; % column # of index.
    axis_st = 1993; % first date for axis of plots 
 end;
 
 
 if data_source == 2; 
    % new data from HFR
    x = load('hfri_historical_yn.txt');
    lbls = {'Distress/Restruct','Merger Arb','Equity Mkt Neutral','Quant Direction','Short',...
          'Emerg Mkts','Equity Hedge','Event-Driven','FOF Conservative','FOF Diversified','FOF Mkt Defensive',...
          'FOF Strategic','FOF Composite','HF Composite','Macro','Macro Syst. Divsfy',...
          'RV Total','RV Fixed & Conv Arb','RV Fixed Inc-Corp','RV Multi-Strat'} ;
    %    
    % lbls = {'Distressed/Restructuring','Merger Arbitrage','Equity Market Neutral','Quantitative Directional','Short Bias',...
    %                  1                       2                   3                           4                 5       
    %        'Emerging Markets','Equity Hedge','Event-Driven','FOF: Conservative','FOF: Diversified','FOF: Market Defensive',...
    %                    6               7            8                  9                   10               11                    
    %        'FOF: Strategic','Fund of Funds Composite','Fund Weighted Composite','Macro (Total)','Macro: Systematic Diversified',...
    %                   12                13                            14             15                      16
    %        'Relative Value (Total)','RV: Fixed Income-Convertible Arbitrage','RV: Fixed Income-Corporate','RV: Multi-Strategy'} ;
    %                   17                                18                              19                      20         
    subset = [14 2 3 4 6 13 15 19] ; % interesting ones for plotting

    THF = size(x,1); 
    r = x(:,3:end) ; % returns; cut off date field.
    plotdate = x(:,2)/12+x(:,1);   
    plotdate_nav = [plotdate(1)-1/12; plotdate]; %??
    allindx = 14; % index(column) of all funds index
    axis_st = 1989; % first date for axis of plots 
 end;
    
 nav = cumprod([100*ones(1,size(r,2)); 1+r/100]); 
 
 % form date vectors for x axis of plots over time

 % **************************
 % Factor data. Option 1, create from underlying files, use when updating
 % data
 % *************************
 
if create_data;
 
 x =  load('F-F_Research_Data_Factors_monthly.txt');
 if data_source == 1; 
    ff_start_date = 199401;
 end;
 if data_source == 2; 
    ff_start_date = 199001;    
 end;
 sti = find(x(:,1) == ff_start_date);
 ff_dates = x(sti:sti+THF-1,1);
 rmrf = x(sti:sti+THF-1,2); % allows for factors longer than returns
 smb = x(sti:sti+THF-1,3);
 hml = x(sti:sti+THF-1,4);
 rf = x(sti:sti+THF-1,5);
 
 x = load('F-F_Momentum_Factor_monthly.txt'); 
 sti = find(x(:,1) == ff_start_date);
 umd = x(sti:sti+THF-1,2); 
  
 
 if 0
     % Data that are no longer in use, preserved in case I want to go
     % back and get them to work.
     % load ibbotson data
     % I haven't been able to update this data, so we have to use finance.yahoo data below instead.
     
     x = load('sbbi.csv');
     x = x(817:end,:);
     rlcorp = x(:,3); % ibbotson long term corporate
     rf_sbbi = x(:,8); % ibbotson 30 day t bill
     cpi = x(:,9);
     rlgovt =  x(:,4); % ibbotson long term govt
     rigovt = x(:,6); % intermediate government
     
     % corp bond returns II
     x= load('VWEHX_2011.csv');
     x = x(end:-1:1,:); % reverse order
     if data_source == 1;
         sti = 55; % feb 1 1994 is return for jan 1994
     end;
     if data_source == 2;
         sti = 1; % needs fixing
     end;
     
     rlcorp = x(sti:end,9)-rf;
     % put write from CBOE
     % this is stil here to look at, but I've abandoned it.
     % though a return, it's at the money puts so ends up looking just like the
     % market (delta) along with a little bit of short-vol. But even there you
     % don't see much different from the market. We really want deep OTM puts,
     % and maybe held to expiration.
     
     x = load('putwrite_new.csv');
     if data_source == 1;
         sti = 92;
     end;
     if data_source == 2;
         sti = 1;
     end;
     
     putwrite = x(sti:end,5)-rf(1:size(x(sti:end,:),1));
     Tput = size(putwrite,1);
     
     
 end;
 
  % govt bond returns
  x = load('bill_inflation.txt');
  %     1            2             3             4             5             6            7              8            9             10            11                                                                                                                            
  %   caldt        b30ret        b20ret        b10ret         b7ret         b5ret         b2ret         b1ret        t90ret        t30ret        cpiret
   if data_source == 1; 
     bill_start_date = 19940131;
   end;
   if data_source == 2; 
     bill_start_date = 19900131;    
   end;
 
  sti = find(x(:,1) == bill_start_date); 
  rlgovt = 100*(x(sti:end,3));% 20 year bond less treasury; 
  rtb    = 100*(x(sti:end,10)); % One month treasury -- wecould use rf, but keep data sources together;   
  rlgovt = rlgovt - rtb; % keep excess return 20 year bond less treasury; 
  
  
  % Attempt to use BAA AAA yields. 
  x = load('BAA_AAA.txt');
  indx = find((x(:,2) == 01)&(x(:,3)==1994));
  x = x(indx-1:end,:);
  AAA = x(:,4);
  BAA = x(:,5);
  rBAA = -7*(BAA(2:end)-BAA(1:end-1) ) + 1/12*(BAA(2:end));
  rAAA = -7*(AAA(2:end)-AAA(1:end-1) ) + 1/12*(AAA(2:end));
 
  rlcorp = rBAA-rAAA; 
  % approximate percent return by change in yield assuming 7 year duration
  
    
    % -------------------------------------
    % output data to file for students
    % ------------------------------------
    
    fid = fopen('..\data\hedge_factors.txt','w');
    fprintf(fid,'%% Factor data for use with hedge fund problem set. Created by ps8\\hedge_index.m \n');
    fprintf(fid,'%% Fama French factors from French website. Term is CRSP 30 year govt bond return - 1 month rate.\n'); 
    fprintf(fid,'%% Def is the change in yield on Fed BAA- AAA\n');
    %fprintf(fid,'%% Def is the vanguard high-yield fund from finance.yahoo.com.\n');
    %fprintf(fid,'%% Put is the Agarwal and Naik put index, less rf, and scaled by -1/20. DO NOT USE for problem set 8\n');
    fprintf(fid,'%% rf is the one month bill rate. All other factors are net of rf, zero cost portfolios.\n');
    fprintf(fid,'%% All returns are monthly percent returns.\n'); 
    fprintf(fid,'%% %10.9s %10.9s %10.9s %10.9s %10.9s %10.9s %10.9s %10.9s \n', 'date','rf', 'rmrf','hml','smb','umd','term','def');
    for t = 1:size(rmrf,1);
        fprintf(fid,' %10.0i  %10.6f %10.6f %10.6f %10.6f %10.6f',ff_dates(t,:),rf(t),rmrf(t),hml(t),smb(t),umd(t));
        if t <= size(rlgovt,1);
            fprintf(fid,' %10.6f',rlgovt(t));
        else
            fprintf(fid,' %10.6s','NaN');
        end;
        if t <= size(rlcorp,1);
            fprintf(fid,' %10.6f',rlcorp(t));
        else
            fprintf(fid,' %10.6s','NaN');
        end;
        fprintf(fid,'\n');
    end
    fclose(fid)
    
end; % ends if create data

% *********************
% Problem set starts
% Load data if already created -- this is what your program will look like
% ********************

if create_data == 0; % load from file
    x = load('hedge_factors.txt');
    ff_dates = x(:,1:3); 
    rf = x(:,2); 
    rmrf = x(:,3);
    hml = x(:,4); 
    smb = x(:,5); 
    umd = x(:,6); 
    rlgovt = x(:,7); 
    rlcorp = x(:,8);
end;

if any(isnan(rlgovt))
    disp('warning nans in rlgovt')
    rlgovt = rlgovt(~isnan(rlgovt));                %      remove NaNs
end

rm = rmrf + rf; 
TFF = size(rm,1);
% form HF excess returns
rx = r - rf*ones(1,size(r,2));
Tterm = size(rlgovt,1); 
Tcorp = size(rlcorp,1); 
    
  
  
 % ------------------
 % Problem set pause. Graphs of 
 % one year returns for my own use. Skip this section 
 % ------------------
 
 if graphs; % comments out this section 
     
 % short sample plot of index to make hf look good
 shortsmp = find(plotdate<2008);
 figure;
 plot(plotdate(shortsmp),100*(exp(filter(ones(12,1),1,log(1+r(shortsmp,allindx)/100)))-1),'-r','linewidth',2);
 hold on
 plot(plotdate,0*plotdate,'-k'); 
 legend('HF',3 );
 axis([axis_st maxdate -40 40]); 
 title('Hedge Fund Index Annual Returns to 2008'); 
 print -depsc2 hedge_return_history_old.eps; 

 % longer sample to see risk 
 figure;
 plot(plotdate,100*(exp(filter(ones(12,1),1,log(1+r(:,allindx)/100)))-1),'-r','linewidth',2);
 hold on
 plot(plotdate,0*plotdate,'-k'); 
 legend('HF',3 );
 axis([axis_st maxdate -40 40]); 
 title('Hedge Fund Index Annual Returns'); 
 print -depsc2 hedge_return_history_a.eps; 

 
 % index with rmrf and def
 figure;
 plot(plotdate,filter(ones(12,1),1,r(:,allindx)),'-r','linewidth',2);
 hold on 
 plot(plotdate,filter(ones(12,1),1,rm(:,1)),'-b','linewidth',2);
 plot(plotdate,filter(ones(12,1),1,rlcorp+rf),'-g','linewidth',2);     
 plot(plotdate,0*plotdate,'-k'); 
 legend('HF','Rm','Def',3 );
 axis([axis_st maxdate -40 40]); 
 title(char(lbls(allindx)')); 
 print -depsc2 hedge_return_history_b.eps; 

 % everyone with rmrf
 for i = 1:size(r,2); 
     figure
     plot(plotdate,100*(exp(filter(ones(12,1),1,log(1+r(:,i)/100)))-1),'-r','linewidth',2);
     hold on 
     plot(plotdate,100*(exp(filter(ones(12,1),1,log(1+rm(:,1)/100)))-1),'-c','linewidth',2);
     plot(plotdate,0*plotdate,'-k'); 
     legend('Hedge Fund','Market',3 );
     axis([axis_st maxdate -50 60]); 
     title(char(lbls(i)')); 
     print('-depsc2',['one_yr_returns_' char(lbls(i)) '.eps' ])
 end;
 
 % one year factor returns
 if create_data == 1; % I didn't save the put data as it was useless. 
     f = {hml smb umd rlgovt rlcorp};
     flbl = {'hml','smb', 'umd', 'term', 'def' };
 end
 if create_data == 0;
     f = {hml smb umd rlgovt rlcorp};
     flbl = {'hml','smb', 'umd', 'term', 'def'};
 end
 
 for i = 1:size(f,2); 
     figure
     f1 = f{i}; 
     T1 = size(f1,1);
     plot(plotdate(1:T1),100*(exp(filter(ones(12,1),1,log(1+(f1+rf(1:T1))/100)))-1),'-r','linewidth',2);
     hold on 
     plot(plotdate,100*(exp(filter(ones(12,1),1,log(1+(rmrf+rf)/100)))-1),'-b','linewidth',2);
     plot(plotdate,0*plotdate,'-k'); 
     legend([char(flbl(i))],'Rm',3 );
     axis([axis_st maxdate -50 60]); 
     print('-depsc2',['one_yr_factor_returns_' char(flbl(i)) '.eps' ])
     %title(char(lbls(i)')); 
 end;
 
 
 end; % ends if graphs to comment out this section. 
 

 % ************************************
 % Problem set starts again 
 % ***********************************
 
 % For regression analysis, chop everyone down to the shortest series 
  
 % --------------------
 % autocorrelations 
 % -------------------- 
 
 disp('serial correlation of hedge fund returns and factors'); 

 % -------------------------
 % simple version limited to one lagged variable -- this is what you do in your problem set
 % -------------------------

 T = size(r,1);
 
 f = [ rmrf hml smb umd rlcorp ...
       [rlgovt; NaN*ones(T-size(rlgovt,1),1) ]]; % in some data updates rlgovt is shorther than others. not now. 
 flbl = {'rmrf','hml','smb', 'umd', 'def','term'}; 
 biglbl = [flbl lbls]; 

 fprintf(' %20s %8.8s %8.8s %8.8s %8.8s \n', 'return','b = rho','t','R2','v ratio')
 bigr = [f rx];
 N2 = size(bigr,2);
 for i = 1:N2;
    rl = bigr(~isnan(bigr(:,i)),i);
    X = [ones(size(rl,1)-1,1) rl(1:end-1)];
    Y = rl(2:end); 
    b = X\Y; 
    e = Y-X*b; 
    s2 = var(e); 
    R2(i) = 1-s2./var(Y); 
    se = (diag(inv(X'*X))*s2).^ 0.5; 
    t = b./se;
    long = filter(ones(12,1),1,log(1+(rl)/100));
    var_long = var(long(13:end));
    var_short = var(log(1+(rl(13:end))/100));
    fprintf(' %20s %8.2f %8.2f %8.3f %8.3f\n', char(biglbl(i)),b(2),t(2),R2(i),var_long/(12*var_short));
 end; 
 
 disp('  '); 
 
 
 % --------------
 % performance attribution regressions
 % --------------
 
 % first, replicate asness et al. In their dataset; afterwards, and full
 % sample
 if 1; 
     if data_source==1; 
         end_date = 129-48;
     end;
     if data_source == 2; 
         end_date = 129;
     end;
     rx1 = rx(1:end_date,allindx); 
     f = rmrf(1:end_date,:); 
     flbls = {'rmrf'}; 
	 rlbls = lbls(allindx); 
     disp('1. Replicate asness -- capm on overall index in short sample');
     fprintf('Date range from  %1i %2i %4i to %1i %2i %4i \n',HFdates(1,:), HFdates(end_date,:)); 
     asness(rx1,f,rf(1:end_date,:),rlbls,flbls,0,3); 

     %disp('2. no lags and up/down for graph'); 
     % commented out scatterplot not useful
	 %asness(rx1,f,rf,rlbls,flbls,0,0);  
     
     
    disp(' '); 
    disp(' ************************' ); 
    disp(' '); 
    f = [rmrf(end_date+1:end,:)];
    flbls = {'rmrf'};
    disp('2. HF Index, 2001-on sample, CAPM only.');
    fprintf('Date range from  %1i %2i %4i to %1i %2i %4i \n',HFdates(end_date+1,:), HFdates(end,:)); 
    asness(rx(end_date+1:end,allindx),f,rf(end_date+1:end,:),lbls(allindx),flbls,0,3);  

    
    disp(' '); 
    disp(' ************************' ); 
    disp(' '); 
    f = [rmrf];
    flbls = {'rmrf'};
    disp('3. HF Index, full sample, CAPM only'); 
    asness(rx(:,allindx),f,rf,lbls(allindx),flbls,0,3);  

    
 end; 
 
%  All styles full sample, CAPM only,

disp(' '); 
disp(' ************************' ); 
disp(' '); 
f = [rmrf];
flbls = {'rmrf'};
disp('3. full sample, CAPM only'); 
asness(rx(:,subset),f,rf,lbls(subset),flbls,0,3);  


% now, same stuff on full sample, include all the style indices
disp(' '); 
disp(' ************************' ); 
disp(' ');  
T1 = min([size(rmrf,1);size(rlgovt,1);size(rlcorp,1)]);
f = [rmrf(1:T1) hml(1:T1) smb(1:T1) umd(1:T1) rlgovt(1:T1) rlcorp(1:T1) ]; % putwrite seems not to be working -- get better data
flbls = {'rmrf','hml','smb', 'umd', 'term','def' };% 'put'
disp('5. full sample, all styles, 3 lags'); 
asness(rx(1:T1,subset),f,rf(1:T1),lbls(subset),flbls,0,3);  


% ------------------------------------------
% Investigate AN puts  
% ------------------------------------------

if (create_data == 1) & 0;
    %  put  sample, CAPM only
    disp(' ');
    disp(' ************************' );
    disp(' ');
    T1 = min(size(rmrf,1),size(putwrite_an,1));
    f = [rmrf(1:T1)];
    flbls = {'rmrf'};
    disp('7. put sample, CAPM');
    asness(rx(1:T1,subset),f,rf(1:T1),lbls(subset),flbls,0,3);
    
    % put sample, puts only
    disp(' ');
    disp(' ************************' );
    disp(' ');
    T1 = min(size(rmrf,1),size(putwrite_an,1));
    f = [putwrite_an(1:T1)];
    flbls = {'put'};
    disp('7. put sample,  AN puts');
    asness(rx(1:T1,subset),f,rf(1:T1),lbls(subset),flbls,0,3);
    
    
    %  put sample, CAPM and puts only
    disp(' ');
    disp(' ************************' );
    disp(' ');
    T1 = min(size(rmrf,1),size(putwrite_an,1));
    f = [rmrf(1:T1) putwrite_an(1:T1)];
    flbls = {'rmrf','put'};
    disp('7. put sample, CAPM and AN puts');
    asness(rx(1:T1,subset),f,rf(1:T1),lbls(subset),flbls,0,3);
    
    % put sample, puts only
    disp(' ');
    disp(' ************************' );
    disp(' ');
    T1 = min(size(rmrf,1),size(putwrite_an,1));
    f = [putwrite(1:T1)];
    flbls = {'put'};
    disp('7. put sample,  CBOE puts');
    asness(rx(1:T1,subset),f,rf(1:T1),lbls(subset),flbls,0,3);
    
    
    %  put sample, CAPM and puts only
    disp(' ');
    disp(' ************************' );
    disp(' ');
    T1 = min(size(rmrf,1),size(putwrite_an,1));
    f = [rmrf(1:T1) putwrite(1:T1)];
    flbls = {'rmrf','put'};
    disp('7. put sample, CAPM and CBOE puts');
    asness(rx(1:T1,subset),f,rf(1:T1),lbls(subset),flbls,0,3);
    
end;
 
stop; 

% *********************
% Ideas not used in problem set
% ********************************

  

% ****************
 % plot distributions -- non-normal? 
 % these actually show some tremendous non-normalities. Fixed income arb, for example
 % seems to have totally sold off the top end (written call options). But I can't find 
 % a sexy way to show this unless you're really good at reading normal distribution graphs
 
 figure; 
 for i = 1:size(r,2); 
    subplot(3,5,i);
    r2 = r(~isnan(r(:,i)),i);
    plot(sort(r2),(1:size(r2,1))/size(r2,1));
    hold on; 
    pts = (min(r2)-3:0.01:max(r2)+3)';
   plot(pts,cdf('normal',pts,mean(r2),std(r2)),'-r')
    title(char(lbls(i))); 
    axis([min(r2)-3 max(r2)+3 0 1 ]); 
 end; 
 
 % --------------
 %normplots
 % --------------
 
 % examples
 
 samples = 150; % number of points in normplot
 figure; 
 % normal
 subplot(2,2,1); 
 x = (-3:0.01:3); 
 plot(x,normpdf(x)); 
 title('normal'); 
 subplot(2,2,2); 
 normplot(randn(samples,1)); 
 title('normal');
 
 % uniform
 subplot(2,2,3);
 x = (-0.5:0.01:1.5);
 plot( x, (x>=0)&(x<1)); 
 axis([-0.2 1.2 -0.2 1.2]); 
 title('uniform'); 
 subplot(2,2,4); 
 normplot(rand(samples,1)); 
 title('uniform'); 
 print -depsc2 normplot_example_1.eps; 
 
 % lognormal
 figure; 
 subplot(2,2,1); 
 x = (0:0.02:8);
 y = lognpdf(x,0,1);
 plot(x,y); 
 title('lognormal'); 
 axis([0 8 -inf inf]); 
 subplot(2,2,2); 
 y = exp(randn(samples,1)); 
 y = y(y<5); % gets rid of occasional huge ones that make graph hard to see
 normplot(y); 
 title('lognormal'); 

 % mixture of normals
 subplot(2,2,3); 
 x = (-5:0.01:5)'; 
 plot(x,[0.5*normpdf(x) 0.5*normpdf(x/2)],'k--');
 hold on; 
 plot(x,0.5*normpdf(x)+ 0.5*normpdf(x/2)); 
 title('normal mixture'); 
 subplot(2,2,4); 
 normplot([randn(samples/2,1); 2*randn(samples/2,1)]); 
 title('normal mixture'); 
 print -depsc2 normplot_example_2.eps; 
  
 % --------------
 % data 
 % -------------
 
 figure; 
 subplot(3,5,1); 
 r2 = rmrf+rf;
 normplot(log(1+r2/100)); 
 title('rm')
 
 for i = 1:size(r,2); 
    subplot(3,5,i+1);
    r2 = r(~isnan(r(:,i)),i);
    normplot(log(1+r2/100))
    title(char(lbls(i))); 
 end; 
 
 print -depsc2 normplot_data.eps; 
 
 
  
 for i = 1:size(r,2); 
  figure;
  plot(plotdate, rmrf,'LineWidth',2); 
  hold on; 
  plot(plotdate,r(:,i)-rf,'-r','LineWidth',2); 
  legend('rmrf',char(lbls(i)),3); 
  axis([1994 2005 -inf inf]); 
  print('-depsc2',['return_' char(lbls(i)) '.eps' ])
 end; 
 

 
 % -------------------------
 % following is the same thing with multiple regressions, many lags of returns. It ended up not making any difference, so not assinged inthe
 % problem set 
 % -------------------------
 
 if 0 ; % use to include or exclude this section from results
 maxlag = 4; % # of lags used to forecast returns.  
 
 if maxlag == 1; 
        fprintf(' %20s %8.2s % 8.2s \n', 'return','b = rho','t')
 end; 
 for i = 1:size(f,2); 
    bv = zeros(maxlag,maxlag);
    tv = bv; 
    R2 = zeros(1,maxlag);
    for j = 1:maxlag;  
        X = f(j:T-1,i); 
        for k = 2:j; 
            X = [X f(j+1-k:T-k,i)]; 
        end; 
        X = [ones(T-j,1) X];
        b = X\f(j+1:T,i); 
        bv(1:j,j) = b(2:end); 
        e = f(j+1:T,i)- X*b;
        se = diag(inv(X'*X)*var(e));
        se = se(2:end).^0.5;  
        tv(1:j,j) = b(2:end)./se; 
        R2(j) = 1-var(e)/var(f(j+1:T,i)); 
    end;
    
    if maxlag == 1; 
        fprintf(' %20s %8.2f %8.2f \n', char(flbl(i)), bv, tv)
    else; 
        fprintf(' %20s \n',char(flbl(i))); 
        fprintf(' r(t-1) ')
        fprintf(' %8.2f ', bv(1,:)); 
        fprintf('\n');
        if maxlag > 1; 
            fprintf(' r(t-2)  %8.2s ',' '); 
            fprintf(' %8.2f  ', bv(2,2:end));
            fprintf(' \n');
        end; 
        if maxlag > 2; 
            fprintf(' r(t-3)  %8.2s  %8.2s ',' ',' '); 
            fprintf(' %8.2f  ', bv(3,3:end));
            fprintf(' \n');
        end; 
        if maxlag > 3; 
            fprintf(' r(t-4) %8.2s %8.2s %8.2s %8.2f \n',' ',' ',' ',bv(4,4)); 
         end; 
    
        fprintf('     R2 ')
        fprintf( ' %8.2f ', R2);
        fprintf('\n')
        fprintf('root R2 ')
        fprintf( ' %8.2f ', R2.^0.5);
        fprintf('\n')
    end; 
 end;    
    

 % autocorrelation of hedge fund excess returns
    
 
for i = 1:size(r,2); 
    bv = zeros(maxlag,maxlag);
    R2 = zeros(1,maxlag);
    for j = 1:maxlag;  
        X = rx(j:T-1,i); 
        for k = 2:j; 
            X = [X rx(j+1-k:T-k,i)]; 
        end; 
        X = [ones(T-j,1) X];
        b = X\rx(j+1:T,i); 
        bv(1:j,j) = b(2:end); 
        e = rx(j+1:T,i)- X*b; 
        se = diag(inv(X'*X)*var(e));
        se = se(2:end).^0.5;  
        tv(1:j,j) = b(2:end)./se; 
       
        R2(j) = 1-var(e)/var(rx(j+1:T,i)); 
    end;
    
    if maxlag == 1; 
        fprintf(' %20s %8.2f % 8.2f \n', char(lbls(i)), bv, tv)
    else; 
        fprintf(' %20s \n',char(lbls(i))); 
        fprintf(' r(t-1) ')
        fprintf(' %8.2f ', bv(1,:)); 
        fprintf('\n');
        if maxlag > 1; 
            fprintf(' r(t-2)  %8.2s ',' '); 
            fprintf(' %8.2f  ', bv(2,2:end));
            fprintf(' \n');
        end; 
        if maxlag > 2; 
            fprintf(' r(t-3)  %8.2s  %8.2s ',' ',' '); 
            fprintf(' %8.2f  ', bv(3,3:end));
            fprintf(' \n');
        end; 
        if maxlag > 3; 
            fprintf(' r(t-4) %8.2s %8.2s %8.2s %8.2f \n',' ',' ',' ',bv(4,4)); 
         end; 
    
        fprintf('     R2 ')
        fprintf( ' %8.2f ', R2);
        fprintf('\n')
        fprintf('root R2 ')
        fprintf( ' %8.2f ', R2.^0.5);
        fprintf('\n')
    end; 
 end;    
 end;   % ends "if 1" that started this whole mess above    
 
 
  % ************************************
  % plots of cumulative returns over time
  % superseded by annual return plots. Put if 1 in if you want to see them.
  % 
  % ************************************
  
  if 0;
      
      Tterm = size(rtb,1);
      rlgovt_nav = 100*cumprod([1;1+(rlgovt+rf)/100]);
      rlcorp_nav = 100*cumprod([1;1+(rf+rlcorp)/100]); % nav of long treasury, (BAA-AAA)
      rput_nav = 100*cumprod([1;1+x(sti:end,5)/100]); % nav of put write strategy
      
      rm_nav = 100*cumprod([1;1+(rmrf+rf)/100]);
      hml_nav = 100*cumprod([1;1+(hml+rf)/100]);
      smb_nav = 100*cumprod([1;1+(smb+rf)/100]);
      umd_nav = 100*cumprod([1;1+(umd+rf)/100]);
      
      
      figure;
      plot(plotdate_nav,rm_nav,'linewidth',2);
      hold on
      plot(plotdate_nav(1:size(rput_nav,1)),rput_nav,'-k','linewidth',2);
      text(plotdate_nav(end),rm_nav(end),'Rm');
      text(plotdate_nav(end),rput_nav(end),'PUT');
      title('CBOE put and market');
      print -depsc2 putwrite_history.eps
      
      % Clearly, we need to do better on this put business!
      
      % ---------------------------------------------------
      % NAV plots -- dropped in favor of 1 year return plots
      % ------------------------------------------------
      
      for i = 1: size(r,2);
          figure; % all
          plot(plotdate_nav,log(nav(:,i)),'linewidth',2);
          hold on;
          text(plotdate_nav(THF),log(nav(THF,i))',lbls(i));
          plot(plotdate_nav,log(rm_nav),'-k','linewidth',2);
          text(plotdate_nav(end),log(rm_nav(end)),'Rm');
          axis([1993.5 2013 -inf inf]);
          print('-depsc2',['NAV_' char(lbls(i)) '.eps' ])
      end;
      
      % NAV factor plots
      % f = [hml smb umd rlgovt rlcorp putwrite];
      f = {hml smb umd rlgovt rlcorp putwrite};
      flbls = {'hml','smb', 'umd', 'term', 'def','put'};
      %fnav = 100*[ones(1,6); cumprod(1+(f+rf*ones(1,6))/100)];
      for i = 1:6;
          fnav = 100*[1; cumprod(1+(f{i}+rf(1:size(f{i},1)))/100)];
          
          figure; % all
          plot(plotdate_nav(1:size(fnav,1)),log(fnav),'linewidth',2);
          hold on;
          text(plotdate_nav(T),log(fnav(end))',flbls(i));
          plot(plotdate_nav,log(rm_nav),'-k','linewidth',2);
          text(plotdate_nav(end),log(rm_nav(end)),'Rm');
          axis([1993.5 2013 -inf inf]);
          print('-depsc2',['F_NAV_' char(flbls(i)) '.eps' ])
      end;
      
      
      % plot them all against the market
      
      for i = 1: size(r,2);
          figure
          plot(plotdate,r(:,i),'-r','linewidth',2);
          hold on
          plot(plotdate,rm(:,1),'-b','linewidth',2);
          plot(plotdate,0*plotdate,'-k');
          legend('HF ','market',3 );
          axis([axis_st 2009 -15 10]);
          title(char(lbls(i)'));
          print('-depsc2',['one_mo_returns_' char(lbls(i)) '.eps' ])
      end;
  end; % ends if to do nav plots or not

  
  
  stop 
  
  % deleted: use agarwal put data to benchmark put writing. Turns out this
  % is not a good data set, focusing on atm puts so mostly just atm implied
  % volatility
  
% ------------------
% Agarwal put data
% ------------------

x = load('Agarwal_optionreturns_0283_0809.csv');
if data_source == 1; 
    x = x(131:end,:);
end;
Topt = size(x,1); % alas usually less than full sampe
option_date = x(:,1:2);   
atm_call = x(:,3)-rf(1:Topt);
otm_call = x(:,4)-rf(1:Topt);
atm_put = x(:,5)-rf(1:Topt);
otm_put = x(:,6)-rf(1:Topt); 

disp('put portfolio excess returns')
fprintf('          %6.6s %6.6s %6.6s \n','mean','std','sharpe');
fprintf(' ATM call %6.2f %6.2f %6.2f \n',mean(atm_call), std(atm_call), mean(atm_call)/std(atm_call))
fprintf(' ATM put  %6.2f %6.2f %6.2f \n',mean(atm_put), std(atm_put), mean(atm_put)/std(atm_put))
fprintf(' OTM call %6.2f %6.2f %6.2f \n',mean(otm_call), std(otm_call), mean(otm_call)/std(otm_call))
fprintf(' OTM put  %6.2f %6.2f %6.2f \n',mean(otm_put), std(otm_put), mean(otm_put)/std(otm_put))

putwrite_an = -otm_put/20; % divisor scales variance to something more reasonable. 
% Figure shows that moneyness difference is too small -- no difference
% between ATM and OTM
Topt = size(putwrite_an,1);

% ---------------------------
% monthly put graphs
% --------------------------

if graphs
    figure;
    plot(1994+1/12:1/12:1994+Topt/12,putwrite_an,'-r','linewidth',2)
    title('AN put');
    print -depsc2 an_put.eps

    figure;
    plot(1994+1/12:1/12:1994+Topt/12,putwrite_an,'-r','linewidth',2)
    hold on
    plot(1994+1/12:1/12:1994+T/12,rmrf,'-b','linewidth',2)
    title('AN put and rmrf')
    print -depsc2 an_put_and_rmrf.eps;
    
    figure
    plot(1994+1/12:1/12:1994+Tput/12,putwrite,'-r','linewidth',2)
    hold on
    plot(1994+1/12:1/12:1994+T/12,rmrf,'-b','linewidth',2)
    title('CBOE put and rmrf');
    print -depsc2 cboe_puts.eps
    
    figure
    plot(1994+1/12:1/12:1994+Topt/12,putwrite_an,'-r','linewidth',2)
    hold on
    plot(1994+1/12:1/12:1994+Tput/12,putwrite,'-b','linewidth',2)
    title('AN put and CBOE put');
    print -depsc2 an_cboe_puts.eps
    
    figure
    plot(1994+1/12:1/12:1994+Topt/12,putwrite_an,'-r','linewidth',2)
    hold on
    plot(1994+1/12:1/12:1994+Topt/12,-atm_put/20,'-b','linewidth',2)
    title('AN OTM put and AN ATM put');
    print -depsc2 an_puts_compared.eps
end
 