%MATLAB SCRIPT FILE s322hw3.m  */
%
%    SOLUTION TO HOMEWORK 3 FOR STATISTICS 322
%    POSTED 9/12/96
%    EXPLORES LAW OF LARGE NUMBERS FOR THE PARETO DISTRIBUTION

iquest = 1 ;     %  1,2,3,4

format compact ;
running = [10 'MATLAB script s322hw3.m, with iquest = ' num2str(iquest) 10]

valpha = [1.5 2 2.5 3 4 10] ;
          %  DEFINE THIS HERE SINCE IT IS COMMON TO ALL PARTS


if iquest == 1 ;    %  THEN DO QUESTION 1

  nobs = 1000 ;
  rand('seed',23398475) ;
  
  malpha = ones(nobs,1) * valpha ;
          %  MATRIX FORM, WITH COMMON COLUMNS

  mpardat = rand(size(malpha)) .^ (1 ./ (1 - malpha)) ;
          %  MATRIX OF PARETO DATA
  
  note = 'Calculate Empirical Probabilities' 
  %  CALCULATE EMPIRICAL PROBABILITIES
  vxcut = [1.1 1.5 2 3 5] ;
  meprob = [] ;
  for icut = 1:length(vxcut) ;
    xcut = vxcut(icut) ;
    meprob = [meprob; mean(mpardat <= xcut)] ;
           %  BUILD UP MATRIX OF EMPIRICAL PROBABILITIES
  end ;
  
  note = 'Calculate Theoretical Probabilities' 
    malpha = ones(length(vxcut),1) * valpha ;
    mxcut = vxcut' * ones(1,length(valpha)) ;
  mtprob = 1 - mxcut .^ (1 - malpha) ;

  cil = meprob - 1.96 * sqrt(mtprob .* (1 - mtprob)) / sqrt(nobs) ;  
  cir = meprob + 1.96 * sqrt(mtprob .* (1 - mtprob)) / sqrt(nobs) ;  
          %  KNOWN p BINOMIAL CONFIDENCE INTERVALS  

  note = '   Output results to file'
  fid = fopen('s322hw3.out','wt') ;
    titlstr1 = 'Output from the MATLAB Script s322hw3.m,   ' ;
    titlstr1 = [titlstr1,date] ;
    titlstr2 = '  Solution to Home Work Assignment 3, ' ;
    titlstr2 = [titlstr2 'due September 12, 1996'] ;
  cntbytes = fprintf(fid,'%1s\n',titlstr1) ;
  cntbytes = fprintf(fid,'%1s\n',titlstr2) ;

    headstr1 = 'Problem 1:  ' ;
  cntbytes = fprintf(fid,'\n\n %1s\n',headstr1) ;

  for ialpha = 1:length(valpha) ;
    alpha = valpha(ialpha) ;
    cntbytes = fprintf(fid,'\n              alpha = %7.3f\n',alpha) ;
      frmtstr = '  xcut    %7.3f    %7.3f    %7.3f    %7.3f    %7.3f\n' ;
    cntbytes = fprintf(fid,frmtstr,vxcut) ;
      frmtstr = '   cil    %7.3f    %7.3f    %7.3f    %7.3f    %7.3f\n' ;
    cntbytes = fprintf(fid,frmtstr,cil(:,ialpha)) ;
      frmtstr = ' tprob    %7.3f    %7.3f    %7.3f    %7.3f    %7.3f\n' ;
    cntbytes = fprintf(fid,frmtstr,mtprob(:,ialpha)) ;
      frmtstr = '   cir    %7.3f    %7.3f    %7.3f    %7.3f    %7.3f\n' ;
    cntbytes = fprintf(fid,frmtstr,cir(:,ialpha)) ;
  end ;

elseif iquest == 2 ;    %  THEN DO QUESTION 2

  vn = [10; 100; 1000;  10000] ;
  rand('seed',98743475) ;
  
  mmean = [] ;
  msd = [] ;
  for in = 1:length(vn) ;
    n = vn(in) ;
    malpha = ones(n,1) * valpha ;
    mpardat = rand(size(malpha)) .^ (1 ./ (1 - malpha)) ;
          %  MATRIX OF PARETO DATA
    mmean = [mmean; mean(mpardat)] ;
    msd = [msd; std(mpardat)] ;
  end ;

  note = '   Output results to file'
  fid = fopen('s322hw3.out','at') ;
    headstr1 = 'Problem 2:   Sample Means' ;
  cntbytes = fprintf(fid,'\n\n\n %1s\n',headstr1) ;
    frmtstr = '  alpha %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f\n' ;
  cntbytes = fprintf(fid,frmtstr,valpha) ;
  cntbytes = fprintf(fid,'%1s\n','    n ') ;
    frmtstr = ' %6.0f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f\n' ;
  cntbytes = fprintf(fid,frmtstr,[vn'; mmean']) ;
    
    headstr1 = 'Problem 2:   Sample Standard Deviations' ;
  cntbytes = fprintf(fid,'\n %1s\n',headstr1) ;
    frmtstr = '  alpha %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f\n' ;
  cntbytes = fprintf(fid,frmtstr,valpha) ;
  cntbytes = fprintf(fid,'%1s\n','    n ') ;
    frmtstr = ' %6.0f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f\n' ;
  cntbytes = fprintf(fid,frmtstr,[vn'; msd']) ;

elseif iquest == 3 ;    %  THEN DO QUESTION 3

  vmean = (valpha - 1) ./ (valpha - 2) ;
  vmean(1) = 0/0 ;   
          %  DELIBERATELY INTRODUCE NAN, SINCE MEAN NOT DEFINED THERE

  vsd = sqrt(((valpha - 1) ./ (valpha - 3)) - vmean.^2) ;
  vsd(1) = 0/0 ;
  vsd(2) = 0/0 ;
  vsd(3) = 0/0 ;

  note = '   Output results to file'
  fid = fopen('s322hw3.out','at') ;
    headstr1 = 'Problem 3:   Theoretical Moments' ;
  cntbytes = fprintf(fid,'\n\n\n %1s\n',headstr1) ;
    frmtstr = '    alpha  %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f\n' ;
  cntbytes = fprintf(fid,frmtstr,valpha) ;
    frmtstr = '     mean  %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f\n' ;
  cntbytes = fprintf(fid,frmtstr,vmean) ;
    frmtstr = 'std. dev.  %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f\n' ;
  cntbytes = fprintf(fid,frmtstr,vsd) ;

elseif iquest == 4 ;    %  THEN DO QUESTION 4

  nset = 500 ;
  nobs = 100 ;
  rand('seed',23889687) ;
  
  vcovp = [] ;
  for ialpha = 3:length(valpha) ;
    alpha = valpha(ialpha) ;
    vt = (alpha - 1) / (alpha - 2) ;
    mpardat = rand(nobs,nset) .^ (1 ./ (1 - alpha)) ;
          %  MATRIX OF PARETO DATA
    vmean = mean(mpardat) ;
    vsd = std(mpardat) ;
    cil = vmean - 1.96 * vsd / sqrt(nobs) ;
    cir = vmean + 1.96 * vsd / sqrt(nobs) ;
          %  ENDPOINTS OF STANDARD C.I.'S
    covp = length(find((cil <= vt)&(vt <= cir))) / nset ;
          %  COVERAGE PROBABILITIES
    vcovp = [vcovp covp] ;
  end ;

  vcpse = sqrt(vcovp .* (1 - vcovp) ./ nset) ;
          %  COVERAGE PROBABILITY STANDARD ERRORS
          %  (CALCULATED ASSUMING BINOMIAL MODEL)

  note = '   Output results to file'
  fid = fopen('s322hw3.out','at') ;
    headstr1 = 'Problem 4:   Coverage Probabilities' ;
  cntbytes = fprintf(fid,'\n\n\n %1s\n',headstr1) ;
    frmtstr = '    alpha   %8.2f %8.2f %8.2f %8.2f\n' ;
  cntbytes = fprintf(fid,frmtstr,valpha(3:6)) ;
    frmtstr = 'Cov. Pro.   %8.4f %8.4f %8.4f %8.4f\n' ;
  cntbytes = fprintf(fid,frmtstr,vcovp) ;
    frmtstr = 'std. err.   %8.4f %8.4f %8.4f %8.4f\n' ;
  cntbytes = fprintf(fid,frmtstr,vcpse) ;


end ;

fclose('all') ;


