%	1 = MSCI Hong Kong	
%	2 = MSCI Singapore	
%	3 = MSCI Brazil	
%	4 = MSCI Argentina	
%	5 = MSCI UK	
%	6 = MSCI Germany	
%	7 = MSCI Canada	
%	8 = MCSI France	
%	9 = MCSI Japan	
%	10 = S&P 500	

%	Monthly data, Jan 1988 to Jan 2002

load countries.dat ;
[n,m] = size(countries) ;
price = countries(:,4:m) ;

log_retn = diff(log(price)) ;
bmu_true = mean(log_retn)' ;
N = length(bmu_true) ;

bOmega_true = cov(log_retn) ;
ibOmega_true = inv(bOmega_true) ;
fsize = 12 ;
ngrid = 1 ;
muP = .012 ;
bone = ones(N,1) ;

A_true = bone'*ibOmega_true*bmu_true ;
B_true = bmu_true'*ibOmega_true*bmu_true ;
C_true = bone'*ibOmega_true*bone ;
D_true = B_true*C_true - A_true^2 ;
bg_true = (B_true*ibOmega_true*bone - A_true*ibOmega_true*bmu_true)/D_true ;
bh_true = (C_true*ibOmega_true*bmu_true - A_true*ibOmega_true*bone)/D_true ;
gg_true = bg_true'*bOmega_true*bg_true ;
hh_true = bh_true'*bOmega_true*bh_true ;
gh_true = bg_true'*bOmega_true*bh_true ;

omegaP_opt =  bg_true + muP*bh_true ;
sigmaP_opt = sqrt(omegaP_opt'*bOmega_true*omegaP_opt) ;

muP_ef = linspace(0,.017,25) ;
omegaP_ef = zeros(N,25) ;

for i = 1:25 ;
omegaP_ef(:,i) =  bg_true + muP_ef(i)*bh_true ;
sigmaP_ef(i) = sqrt(omegaP_ef(:,i)'*bOmega_true*omegaP_ef(:,i)) ;
end ;


rand('seed',55393775) ;
for iter = 1:400 ;

select = ceil((n-1)*rand(n-1,1)) ;
log_retn_boot = log_retn(select,:) ;

bmu_boot = mean(log_retn_boot)' ;
bOmega_boot = cov(log_retn_boot) ;




%	Compute minimum expected return and minimum return SD
ibOmega_boot = inv(bOmega_boot) ;	% Invert Omega - "i" means inverse

A_boot = bone'*ibOmega_boot*bmu_boot ;
B_boot = bmu_boot'*ibOmega_boot*bmu_boot ;
C_boot = bone'*ibOmega_boot*bone ;
D_boot = B_boot*C_boot - A_boot^2 ;
bg_boot = (B_boot*ibOmega_boot*bone - A_boot*ibOmega_boot*bmu_boot)/D_boot ;
bh_boot = (C_boot*ibOmega_boot*bmu_boot - A_boot*ibOmega_boot*bone)/D_boot ;
gg_boot = bg_boot'*bOmega_boot*bg_boot ;
hh_boot = bh_boot'*bOmega_boot*bh_boot ;
gh_boot = bg_boot'*bOmega_boot*bh_boot ;




omegaP_boot =  bg_boot + muP*bh_boot ;
sigmaP_boot = sqrt(omegaP_boot'*bOmega_true*omegaP_boot) ;
muP_boot = bmu_true'*omegaP_boot ;



sigmaP_boot_stor(iter) = sigmaP_boot ;
muP_boot_stor(iter) = muP_boot ;

end ;
figure(1)
clf
subplot(3,1,1)
fsize = 14 ;
p = plot(sigmaP_boot_stor,muP_boot_stor,'.',sigmaP_opt,muP,'o', ...
	sigmaP_ef,muP_ef) ;
set(p(2),'color','black') ;
set(p(2),'markersize',60)
set(p(2),'marker','.') ;
set(p(3),'color','magenta') ;
set(p(3),'linewidth',5) ;
set(gca,'fontsize',fsize) ;
set(gca,'ylim',[.008 .016]) ;
xlabel('standard deviation of return (\sigma_P)','fontsize',fsize) ;
ylabel('expected return (\mu_P)','fontsize',fsize) ;
legend('achieved','optimal','eff. frontier',2)

subplot(3,1,2) ;
hist(sigmaP_boot_stor ./ sigmaP_opt,14) ;
colormap('gray')
xlabel('\sigma_P/\sigma_{P,opt}','fontsize',fsize) ;
ylabel('frequency','fontsize',fsize) ;
set(gca,'fontsize',fsize) ;

subplot(3,1,3) ;
hist(muP_boot_stor ./ .012,14) ;
xlabel('\mu_P/.012','fontsize',fsize) ;
ylabel('frequency','fontsize',fsize) ;
set(gca,'fontsize',fsize) ;


orient tall ;

print count_all03.ps -deps ;