%% load historical price data
try
    price; % do not re-load data if it is already in the workspace
catch
    [price last]=yahoo_prices({'^DJI'...
        'AA' 'AXP' 'BA' 'BAC' 'CAT' 'CSCO' 'CVX' 'DD' 'DIS' 'GE'...
        'HD' 'HPQ' 'IBM' 'INTC' 'JNJ' 'JPM' 'KFT' 'KO' 'MCD' 'MMM'...
        'MRK' 'MSFT' 'PFE' 'PG' 'T' 'TRV' 'UTX' 'VZ' 'WMT' 'XOM'},...
        '31-Jan-2010','31-Jan-2012');
    tickers=gettimeseriesnames(price);
    index_ticker=tickers{1};tickers(1)=[];
% remove holidays and fill in missing values based on the index
    price=resample(price,price.Time(~isnan(price.(index_ticker).Data)));
    dates=getabstime(price);
end
%% fit GARCH, forecast variance, and standardize residuals
obj=@(h,epsilon)sum(log(2*pi*h)+epsilon.^2./h); % quasi-MLE objective
condvar=tscollection(dates,'Name','NGARCH conditional variance');
condvar.TimeInfo=price.TimeInfo;
z=[];p=[];muX=[];stdX=[];
for ticker=tickers;j=ticker{:};
    y=log(price.(j).Data(2:end,:)...
        ./price.(j).Data(1:end-1,:)); % log total returns
    m=zeros(size(y)); % conditional means (zero)
    epsilon=y-m; % residuals
    [params val flag.(j)]=fminsearch(@(params)...
        obj(NGARCH2(epsilon,params),epsilon),[log(.1) log(.8) 0]);
    logL.(j)=-val; % log likelihood
    [h omega.(j) h1 h0]=NGARCH2(epsilon,params);
    condvar=addts(condvar,timeseries([h0;h],dates,'Name',j));
    alpha.(j)=exp(params(1));
    beta.(j)=exp(params(2));
    gamma.(j)=params(3);
    phi=alpha.(j)*(1+gamma.(j)^2)+beta.(j);
    varX=(1-phi^10)/(1-phi)*(h1-omega.(j)/(1-phi))...
        +10*omega.(j)/(1-phi); % ten-day log-return variance forecast
    stdP.(j)=last.(j)*sqrt(varX); % price standard deviation by delta rule
    z=[z epsilon./sqrt(h)]; % standardized residuals
    p=[p;last.(j)]; % array of close prices
    muX=[muX;10*mean(y)]; % array of return sample mean (really naive!)
    stdX=[stdX;sqrt(varX)]; % array of forecasted return std dev
end
clear ticker j y epsilon params val h h0 h1 phi varX
%% estimate concordances of standardized residuals
count=0;
tau=zeros(size(z,2));
for j=1:(size(z,1)-1)
    for k=(j+1):size(z,1)
        conc=sign(z(j,:)-z(k,:));
        tau=tau+conc'*conc;
        count=count+1;
    end
end
tau=tau/count;
tau=tau-diag(diag(tau))+eye(size(tau)); % adjust diagonal for ties
clear count j k conc
%% identify with Kendall tau for a Gaussian copula
rho=sin(pi/2*tau);
%% forecast moments for log-normal model for 'P'
covX=stdX*stdX'.*rho;
EP=p.*exp(muX+diag(covX)/2);
covP=EP*EP'.*(exp(covX)-1);
%% calculate moments for market vector 'M = P - p'
covM=covP;
EM=EP-p;
%% evaluate efficient frontier
c=1E7; % initial investment
alphaSR=c*inv(covM)*EM/(p'*inv(covM)*EM);
%% store the results
for k=1:length(tickers);j=tickers{k};
    sharesSR.(j)=alphaSR(k);
    dollarsSR.(j)=sharesSR.(j)*last.(j);
    weightSR.(j)=dollarsSR.(j)/c;
end
clear j k
sprintf('leverage %f',sum(structfun(@abs,weightSR)))