% This programs was wrotten by Shengying Li, 
%                              ph.d. student of Stony Brook University. 
%                              lshengyi@cs.sunysb.edu
% Date of release: May 13, 2004 
%
%    
% Copyright (c) 2004                                                  by Shengying Li
%
% Permission is granted to use and distribute the program for research or educational usage. 
% You can copy or modificate the program as you want, once this permission notice is 
% preserved on all copies. 
%
% Please read the readme file before you use it.     
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Read input files
% load price of 10 stocks
load stock.txt
% load price of indices 
load index.txt

%initialize variables
Srow0 = size(stock10,1);
Srow = Srow0-1;
Scol = size(stock10,2);
Irow = Srow+1;             
Icol = size(index,2)+1;

% compute return value to stock and index
stockreturn = zeros(Srow, Scol);
indexreturn = ones(Irow, Icol);

% return for stock
for i=1:Srow
    for j=1:Scol
	stockreturn(i,j) = stock10(i,j)/stock10(i+1,j)-1;
     end
end

% Add the title to indexreturn 
indexreturn(1,1) = 0;
for j=2:Icol
	indexreturn(1,j)=j-1;
end

% return for index
for i=2:Irow
     for j=2:Icol
 	indexreturn(i,j) = index(i-1,j-1)/index(i,j-1)-1;
     end
end

%Compute size column by small-large
for i=2:Irow
	indexreturn(i,Icol-1)=indexreturn(i,Icol-1)-indexreturn(i,Icol);
end
Icol=Icol-1;

% Initialize the out matrix
Out=zeros(Scol,Icol+1);

for Si=1:Scol
	Y = stockreturn(:,Si); 
	Aold = indexreturn(:,1:Icol);

	Anew = Aold(2:Irow,:);   
	%initial model
	X = Anew\Y;


	%Calculate the R-square
	Ypredict=Anew*X;
	Ydev=Y-Ypredict;
	R2=1-(std(Ydev)/std(Y))^2;
	R2adj=1-(1-R2)*(Srow-1)/(Srow-Icol-1);

	OldR2adj = R2adj;
	Xold = X;
	Inewcol = Icol;

	flag = 1;

    while flag == 1,
	MaxR2adj=OldR2adj;	

	%Very carefully delete one column;
	Inewcol = Inewcol-1;
	Aold2 = ones(Irow, Inewcol);
	Aold2(1,1) = 0;

	Imax = 0;
	for Idelete = 2:Inewcol+1
		if Idelete > 2
			for j=2:Idelete-1
				Aold2(:,j)=Aold(:,j);		
			end
		end
		if Idelete <Inewcol+1
			for j=Idelete+1:Inewcol+1
				Aold2(:,j-1)=Aold(:,j);
			end
		end
		Anew2=Aold2(2:Irow,:);
		X = Anew2\Y;

		%Calculate the R-square
		Ypredict=Anew2*X;
		Ydev=Y-Ypredict;
		R2=1-(std(Ydev)/std(Y))^2;
		R2adj=1-(1-R2)*(Srow-1)/(Srow-Inewcol-1);

		if R2adj>MaxR2adj
		    	MaxR2adj = R2adj;
  			Imax = Idelete;
		end
	end

	if MaxR2adj > OldR2adj
		% Update the R2adj
		OldR2adj = MaxR2adj;

		% Update the matrix
		if Imax > 2
			for j=2:Imax-1
				Aold2(:,j)=Aold(:,j);		
			end
		end
		if Imax < Inewcol+1
			for j=Imax+1:Inewcol+1
				Aold2(:,j-1)=Aold(:,j);
			end
		end
		Aold = Aold2;
		Xold = X;
	else
		flag = 0;
	end
    end

    %Calculate the covariance matrix
    Anew=Aold(2:Irow,:);
    Ypredict = Anew*Xold;
    Ydev = Y-Ypredict;
    Var = inv(Anew'*Anew)*std(Ydev)^2;
    finalsize=size(Xold,1);
    Aold2= Aold;	
    Iuse = 1;

    for I=2:finalsize
	if( ( Xold(I)/(Var(I,I)^0.5) < -1.96 ) | ( Xold(I)/(Var(I,I)^0.5) > 1.96 ) )
		Iuse = Iuse+1;
		Aold2(:,Iuse)=Aold(:,I);
         end
    end

    Aold3 = Aold2(:,1:Iuse);
    Anew3 = Aold3(2:Irow,:);
    Xnew = Anew3\Y;
    Ypredict = Anew3*Xnew;
    Ydev = Y- Ypredict;
    R2=1-(std(Ydev)/std(Y))^2
    R2adj=1-(1-R2)*(Srow-1)/(Srow-Iuse-1);

% Fill the out matrix
    Remain = Iuse;
    for k=1:Remain
	Out(Si, Aold3(1,k)+1)=Xnew(k);
    end
    Out(Si,Icol+1) = R2adj;

end

% Output file:
save riskmodel Out -ascii