%author: ankur kumar
%Principal component analysis I - finding out if there is a relation between
%the CHANGE in "closing price" of several stocks and thus reducing down to the stocks
%whose pattern are closely related. Hopefully, create a movie file.
%pca source: http://csnet.otago.ac.nz/cosc453/student_tutorials/principal_components.pdf
%load all the data - the closing price
%format: date open high low close volume

%load all the closing price data - the change in closing price
%October 29 2006 - October 28 2007 : 250 points
%parallel coordinates graphs are for the derivatives of the closing prices
%at a given time t

%USD-YEN when this goes down, highlight all the stocks that are down for
%this time period, and when this goes up highlight all the stocks up

function main()
%load the data and calculate derivative, using central point diff formula
load 'C:\Documents and Settings\labuser\Desktop\data\aapl.csv';
aapl = ddt(flipud(aapl(:,5)));
load 'C:\Documents and Settings\labuser\Desktop\data\amzn.csv';
amzn = ddt(flipud(amzn(:,5)));
load 'C:\Documents and Settings\labuser\Desktop\data\ba.csv';
ba = ddt(flipud(ba(:,5)));
load 'C:\Documents and Settings\labuser\Desktop\data\bidu.csv';
bidu = ddt(flipud(bidu(:,5)));
load 'C:\Documents and Settings\labuser\Desktop\data\bsc.csv';
bsc = ddt(flipud(bsc(:,5)));
%dow jones needs to be scaled for better viewing
load 'C:\Documents and Settings\labuser\Desktop\data\dji.csv';
dji = ddt(flipud(dji(:,5)));
% dji = dji./(max(dji)-min(dji));%sclae dow jones because of high values
dji = dji./10^2; %scale dow jones because of big numbers
load 'C:\Documents and Settings\labuser\Desktop\data\goog.csv';
goog = ddt(flipud(goog(:,5)));
load 'C:\Documents and Settings\labuser\Desktop\data\grmn.csv';
grmn = ddt(flipud(grmn(:,5)));
load 'C:\Documents and Settings\labuser\Desktop\data\gs.csv';
gs = ddt(flipud(gs(:,5)));
load 'C:\Documents and Settings\labuser\Desktop\data\isrg.csv';
isrg = ddt(flipud(isrg(:,5)));
load 'C:\Documents and Settings\labuser\Desktop\data\leh.csv';
leh = ddt(flipud(leh(:,5)));
load 'C:\Documents and Settings\labuser\Desktop\data\lmt.csv';
lmt = ddt(flipud(lmt(:,5)));
load 'C:\Documents and Settings\labuser\Desktop\data\mer.csv';
mer = ddt(flipud(mer(:,5)));
load 'C:\Documents and Settings\labuser\Desktop\data\nok.csv';
nok = ddt(flipud(nok(:,5)));
load 'C:\Documents and Settings\labuser\Desktop\data\ptr.csv';
ptr = ddt(flipud(ptr(:,5)));
load 'C:\Documents and Settings\labuser\Desktop\data\rimm.csv';
rimm = ddt(flipud(rimm(:,5)));
load 'C:\Documents and Settings\labuser\Desktop\data\slb.csv';
slb = ddt(flipud(slb(:,5)));
load 'C:\Documents and Settings\labuser\Desktop\data\usdyen.csv';
usdyen = ddt(flipud(usdyen(:,5)));
load 'C:\Documents and Settings\labuser\Desktop\data\vlo.csv';
vlo = ddt(flipud(vlo(:,5)));
load 'C:\Documents and Settings\labuser\Desktop\data\xom.csv';
xom = ddt(flipud(xom(:,5)));
%end of loading closing price data.. about 5 sectors
%the entire financial dataset spanning 5 sectors
data = [aapl amzn ba bidu bsc dji goog grmn gs isrg leh lmt mer nok ptr rimm slb usdyen vlo xom];
%the labels for the dataset, these don't change ever
labels = [{'aapl'} {'amzn'} {'ba'} {'bidu'} {'bsc'} {'dji'} {'goog'} {'grmn'} {'gs'} {'isrg'} {'leh'} {'lmt'} {'mer'} {'nok'} {'ptr'} {'rimm'} {'slb'} {'usdyen'} {'vlo'} {'xom'}];
kclust = 5; %for five sectors of finance
map = jet(length(data));
fig=graph(data,1,map,labels,'original data'); %graph original data
ndim = length(data(1,:));
%adjust the data for mean
means = mean(data);
size(data)
adjdata = data-repmat(mean(data),length(data),1);
fig=graph(adjdata,2,map,labels,'mean adjusted data'); %graph mean-adjusted data
size(map)
%calculate the covariances 
covr = cov(adjdata);
[v,d] = eig(covr); %calculate eigenvalues, principle components
%list of sorted variances, clustered variances are related
%variances are diag(covariances)
[list,in] = sort(diag(covr))
idx = kmeans(list,kclust); %number of sectors looked at - clusters
v=v(:,in); %reformat the eigvector matrix - the principle components
            %with the lowest variances in the beginning
            
nndata = (v'*adjdata')'; %final data mean adjusted
size(means)
size(nndata)
nndata = nndata + repmat(means,length(nndata(:,1)),1);% add the mean
fig=graph(nndata,3,map,labels,'final pca data'); %graph final data
size(v)
size(nndata)
fig=graph((v*nndata')',4,map,labels,'full reconstruction'); %full reconstruction no dimension reduction, 
                                % but dimensions are moved
%using covr reduce dimensions from 20 to 10
%get the principle (adjdata) data back with dimension reduction
%number of principal components to use, sigeig=1(all)
sigeig=5; 
vv=v([sigeig:length(in)],[sigeig:length(in)]); %take the top sigeig# feature eigenvectors-significant
size(vv)
size(nndata)
rndata = (vv*nndata(:,[sigeig:length(in)])')'; %take final data, the features, to get adjdata and graph it
size(rndata)
fig=graph(rndata,5,map,labels,'dimension reduction reconstruction'); %matches with nndata overall
%write csv files of the matrices
dlmwrite('nndata.csv',nndata);

%group all the related clustered data
rfdata = data(:,in); %reformat the data  according to order
old = idx(1);
numchanges = 1;
yaw = 1;
%dimension reductions
for i=1:length(idx)
    %everytime there is a change of cluster number, add those vectors to
    %separate matrix
    if(idx(i)~=old)
        %rearrange data into clusters with appropriate labels
        clusters.data{numchanges} = rfdata(:,[yaw:i-1]');
        clusters.labels{numchanges}= [labels(in(yaw:i-1))];
        old = idx(i);
        yaw = i;
        numchanges = numchanges + 1
    end
end
%list cluster and its labels
clusters.data{numchanges}=rfdata(:,[yaw:i]);
clusters.labels{numchanges}=[labels(in(yaw:i))];
graph(clusters.data{1},6,map,clusters.labels{1},'cluster 1')
graph(clusters.data{2},7,map,clusters.labels{2},'cluster 2')
graph(clusters.data{3},8,map,clusters.labels{3},'cluster 3')
graph(clusters.data{4},9,map,clusters.labels{4},'cluster 4')
graph(clusters.data{5},10,map,clusters.labels{5},'cluster 5')

%write all clusters to file for grandtour
dlmwrite('k1.csv',clusters.data{1});
dlmwrite('k2.csv',clusters.data{2});
dlmwrite('k3.csv',clusters.data{3});
dlmwrite('k4.csv',clusters.data{4});
dlmwrite('k5.csv',clusters.data{5});
%effects of usd-yen, dow jones, fft
fig=graphagainst(data,11,18,labels,'usd-yen effect'); %usd yen
fig=graphagainst(data,12,6,labels,'dow jones effect'); %dji
fig = graphfft(data,13,labels,'magnitude of n-dim FFT');


%calculate derivative f'(n) = (f(n+1) - f(n-1))/2
function tarr=ddt(arr)
tarr = []; %temporary array
for i=2:length(arr)-1
    tarr=[tarr; (arr(i+1)-arr(i-1))/2];
end

%graph in parallel coordinates with arr(data),figure(num), colormap(m),
%title (t)
function fig=graph(arr,num,m,labs,t)
fig=figure(num);
nd = length(arr(1,:));
hold on
for i=1:length(arr)
    plot(arr(i,:),'color',m(i,:));
end
title(t);
%the axes for parallel coords
axisx = [1;nd];
axisy = [0;0];

label = [];
for i = 1:nd
    axisx = [axisx,[i;i]];
    axisy = [axisy,[-35;35]];
    if(isempty(labs))
        label = [label,{sprintf('x%d',i)}];
    else
        label = labs;
    end
end
line(axisx, axisy, 'color', 'b');           %draw the axes
set(gca, 'xtick',[1:1:nd]);
set(gca, 'xticklabel', label);
hold off

%graph in parallel coordinates with arr(data),figure(num), colormap(m)
%graph the rates of change based on change on a data category (18th in data array),
%dji (6 in data array), ag is against variable
function fig=graphagainst(arr,num,ag,labs,t)
fig=figure(num);
nd = length(arr(1,:));
m = jet(8); %six levels of change can be displayed
dat = arr(:,ag);
hold on
for i=1:length(arr)
      if(dat(i)>0)
          plot(arr(i,:),'color',m(1,:));
      else
          plot(arr(i,:),'color',m(5,:));
      end
end
title(t);
%the axes for parallel coords
axisx = [1;nd];
axisy = [0;0];

label = [];
for i = 1:nd
    axisx = [axisx,[i;i]];
    axisy = [axisy,[-35;35]];
    if(isempty(labs))
        label = [label,{sprintf('x%d',i)}];
    else
        label = labs;
    end
end
line(axisx, axisy, 'color', 'b');           %draw the axes
set(gca, 'xtick',[1:1:nd]);
set(gca, 'xticklabel', label);
hold off

%graph in parallel coordinates with arr(data),figure(num), colormap(m)
%graph the fft in n-dimesions, high frequencies beat the market
function fig=graphfft(arr,num,labs,t)
fig=figure(num);
nd = length(arr(1,:));
kclust = 6;
m = jet(kclust); %six levels of change can be displayed
datfft = [];
for i=1:nd
    datfft = [datfft,abs(fft(arr(:,i)))];
end
idx = kmeans(datfft,kclust); %for coloring purposes
hold on
for i=1:length(datfft)
    plot(datfft(i,:),'color',m(idx(i),:));
end
title(t);
%the axes for parallel coords
axisx = [1;nd];
axisy = [0;0];

label = [];
for i = 1:nd
    axisx = [axisx,[i;i]];
    axisy = [axisy,[-35;35]];
    if(isempty(labs))
        label = [label,{sprintf('x%d',i)}];
    else
        label = labs;
    end
end
line(axisx, axisy, 'color', 'b');           %draw the axes
set(gca, 'xtick',[1:1:nd]);
set(gca, 'xticklabel', label);
hold off
