function hw2q2rate
clc % clear variables, but must come before globals, 
    %    else clears globals too.
fprintf('HW2Q2 code: Convergence of Stock Price SDE solutions');
m = 6;
N = (10.^(2:m+1))';%[100,1000,10000,100000]'; 
fprintf('\nsizeN = [%i,%i]; a size larger than in Ques#2;',size(N));
fprintf('\nN = [%i',N(1,1));
for i=2:m,fprintf(',%i',N(i,1));end,fprintf('];');
T = 2; 
fprintf('\nT = %5.3f;',T);
dt = T./N; 
fprintf('\nsizedt = [%i,%i];',size(dt));
fprintf('\ndt = [%9.3e',dt(1,1));
for i=2:m,fprintf(',%9.3e',dt(i,1));end,fprintf('];');
dtm = dt(m,1);
fprintf('\n\n1: Compute the Reference Maximum Sample Results, i.e., i = m.');
Nm = N(m,1);
tm = 0:dtm:T; 
fprintf('\n\nsizetm = [%i,%i];',size(tm));
mum = zeros(1,Nm+1);
sigmam = zeros(1,Nm+1);
sqrtdt = sqrt(dt);
sqrtdtm = sqrtdt(m,1);
randn('state',5);
Z = randn(1,Nm); % Get Reference Normal Distribution.
S0 = 100;
SSDEm = zeros(1,Nm+1); SSDEm(1,1) = S0;
SEXPm = zeros(1,Nm+1); SEXPm(1,1) = S0;
for j = 1:Nm
    mum(1,j) = 0.23*(1+tm(1,j)/10); % Used for All Values of m.
    sigmam(1,j) = 0.36*(1+tm(1,j)/10); % Used for All Values of m.
    SSDEm(1,j+1) = SSDEm(1,j)*(1+mum(1,j)*dtm+sigmam(1,j)*sqrtdtm*Z(1,j));
    SEXPm(1,j+1) = SEXPm(1,j)*exp((mum(1,j)-sigmam(1,j)^2/2)*dtm ...
        +sigmam(1,j)*sqrtdtm*Z(1,j));
end 
fprintf('\nm = %i; Nm = %i = %7.1e; only Test:',m,Nm,Nm);
fprintf('\nstd(SSDEm-SEXPm) = %9.3e;',std(SSDEm-SEXPm));
%
fprintf('\n\n2: Compute Standard Deviation Results:');
fprintf('\n\nStandard Deviation Convergence Test:');
fprintf('\n \t \t     a)\t\t     b.1)\t        b.2)');
fprintf('\n  i\t N\t std(SSDE-SEXP)\t std(SSDE-SEXACT)   std(SEXP-SEXACT)');
LogStda = zeros(m,1); LogStdb1 = zeros(m,1);
for i = 1:m % just a check at m
    SSDE = zeros(1,N(i,1)+1); SSDE(1,1) = S0;
    SEXP = zeros(1,N(i,1)+1); SEXP(1,1) = S0;
    SEXACT = zeros(1,N(i,1)+1); SEXACT(1,1) = S0;
    for j = 1:N(i,1)
        jm = (j-1)*Nm/N(i,1)+1; % Get Reference Data at Same Times,
        %  similar to sampling one market for different time subsets each sample! 
        SSDE(1,j+1) = SSDE(1,j)*(1+mum(1,jm)*dt(i,1) ...
            +sigmam(1,jm)*sqrtdt(i,1)*Z(1,jm));
        SEXP(1,j+1) = SEXP(1,j)*exp((mum(1,jm)-sigmam(1,jm)^2/2)*dt(i,1) ...
            +sigmam(1,jm)*sqrtdt(i,1)*Z(1,jm));
        SEXACT(1,j+1) = SEXPm(1,j*N(m,1)/N(i,1)+1);
    end
    Stda = std(SSDE-SEXP); LogStda(i,1) = log(Stda);
    Stdb1 = std(SSDE-SEXACT); LogStdb1(i,1) = log(Stdb1);
    fprintf('\n  %i\t%7.1e\t %9.3e\t %9.3e\t    %9.3e'...
        ,i,N(i,1),Stda,Stdb1,std(SEXP-SEXACT));  
end
%
fprintf('\n\n3:  Extra Computation of the Convergence Rates.');
fprintf('\n\n Assuming asymptotically, std(Delta[S]) ~ C*(dt(i,1))^beta;');
fprintf('\n   or log(std(Delta[S])) ~ -beta*log(N(1,i)/T)+log(C).');
betaa = zeros(m,1); betab1 = zeros(m,1);
%
fprintf('\n\nLocal Convergence Rates by Linear Interpolation:');
fprintf('\n \t \t   a)\t\t     b.1)');
fprintf('\n  i\t N\t Beta(i)\t Beta(i)');
for i = 2:m
    betaa(i,1) = - (LogStda(i,1)-LogStda(i-1,1))/log(N(i,1)/N(i-1,1));
    betab1(i,1) = - (LogStdb1(i,1)-LogStdb1(i-1,1))/log(N(i,1)/N(i-1,1));
    fprintf('\n  %i\t%7.1e\t %9.3e\t %9.3e'...
        ,i,N(i,1),betaa(i,1),betab1(i,1));
end
fprintf('\nNote: Only Part a) convergence is Useful.');
%
x = -log(N); y = LogStda; % Must be same size.
CovLogNStd = mean((x-mean(x)).*(y-mean(y))); 
% Caution: MATLAB Cov for Matrices.
VarLogN = mean((x-mean(x)).^2);
BetaLS = CovLogNStd/VarLogN; 
% Fit beta = Cov[-log(N),std(SSDe-SEXP)]/Var[-log(N)].
fprintf('\n\nConvergence Rate Fit by Ordinary Least Squares for Part a) Results:');
fprintf('\nCov[-log(N),log(Stda)] = %9.3e; Var[-log(N)] = %9.3e;'...
        ,CovLogNStd,VarLogN);
fprintf('\nConv. Rate Power: beta_a = Cov/Var = %9.3e; Approx. %3.1f;' ...
    ,BetaLS,BetaLS);
%
fprintf('\nNote: Approximately, convergence is O(1/sqrt(N(i,1))),');
fprintf('\n      which is very appropriate for a Normal Distribution');
fprintf('\n      and simple Monte Carlo simulations.');
%
fprintf('\n');
%
% END HW2Q2.m