% this script demonstrates the calibration of a generalized error
% distribution to historical equity index returns
% J. Dodson, 20070918
%% get data from Yahoo! Finance
csv=urlread('http://ichart.finance.yahoo.com/table.csv?s=%5EDJI&a=00&b=5&c=1998&d=08&e=14&f=2007&g=w');
% Dow Jones Industrial Average has symbol ^DJI with caret char %5E
% S&P 500 index is ^GSPC, so use 's=%5EGSPC' in CGI script arglist
% other arguments are start date, end date, and frequency (d,w,m)
% note that the months count 00-11 rather than the usual 01-12
columns=textscan(csv,'%s%*f%*f%*f%*f%*d%f','Delimiter',',','HeaderLines',1);
% the format is yyyy-mm-dd,Open,High,Low,Close,Volume,Adj Close
% use '*' in format string to skip fields
dates=columns{1};
closes=columns{2};
% plot(timeseries(closes,datenum(dates,'yyyy-mm-dd'),'IsDatenum',true,'Name','^DJI'));
%% calculate returns
returns=closes(1:end-1)./closes(2:end)-1; % simple returns
returns_location=mean(returns);
returns_scale=std(returns); % consider other measures of location and scale
standard_returns=(returns-returns_location)/returns_scale;
% hist(standard_returns,100);
%% fit density - see 'GED_density.m' for definition
log_likelihood=@(xi)sum(log(GED_density(standard_returns,xi)));
% fplot(log_likelihood,[0.8 1.6]);
xi_MLE=fminsearch(@(xi)-log_likelihood(xi),1.0); % we have to search for the MIN of the NEGATIVE
returns_density=@(x)GED_density((x-returns_location)/returns_scale,xi_MLE)/returns_scale; % affine transformation
returns_distribution=@(x)quad(returns_density,-1.,x); % CDF
%% calculate Value-at-Risk
VaR_95=-fzero(@(x)returns_distribution(x)-0.05,0.);
VaR_99=-fzero(@(x)returns_distribution(x)-0.01,0.);
%% create random paths
increments=reshape(returns_location+returns_scale*GED_rand(10*52*1000,xi_MLE),10*52,1000);
paths=closes(1)*cumprod(1+increments);
% plot([repmat(closes(52:-1:1),1,20);paths(:,1:20)]);
%% exercise
% what is the probability of the index falling below 10,000?