% White test for heterskedasticity and robust covariance estimates

%Step 1: Load in the data and arrange it in a useful way.

  dta = load('Ex_05.asc');
  y = dta(:,1)/1000;
  n = length(y);
  x = [ones(n,1) dta(:,4)/1000 dta(:,5)/1000 dta(:,3)];
  [n,k] = size(x);

%Step 2: Calculate OLS estimates.

  b = inv(x'*x)*(x'*y);

%Step 3: Calculate summary statistics.

  df = n - k;
  e = y - x*b;
  sse = e'*e;
  s2 = sse/df;
  ybar = sum(y)/n;
  r2 = 1 - (e'*e)/((y - ybar)'*(y - ybar));
  f = (((y - ybar)'*(y - ybar) - sse)/(k-1))/s2;

%Step 3b: Calculate White test.

  e = y - x*b;
  y2 = e.*e;
  x2 = kron(x,ones(1,k)).*kron(ones(1,k),x);
  [u,v,s] = svd(diag(reshape(tril(ones(k,k)),k*k,1)));
  x2 = x2*u(:,1:k*(k+1)/2);
  b2 = inv(x2'*x2)*(x2'*y2);
  e2 = y2 - x2*b2;
  ybar2 = sum(y2)/n;
  w = 1 - (e2'*e2)/((y2 - ybar2)'*(y2 - ybar2));
  w = w*n;

%Step 4: Calculate covariance matrix, standard errors, and t-ratios.

%  c = s2*inv(x'*x);
  c = inv(x'*x)*(x'*((y2*ones(1,k)).*x))*inv(x'*x); % White cov. estimator
  se = sqrt(diag(c));
  tr = b./se;
  pv = 2*(1 - cdf('t',abs(tr),df));

%Step 5: Print results.

  fprintf('\rdf = %d\r', df);
  fprintf('sse = %14.5f\r', sse);
  fprintf('s2 =  %14.5f\r', s2);
  fprintf('r2 =  %14.5f\r', r2);
  fprintf('f =   %14.5f, p-value = %14.5f\r', [f (1 - cdf('f',f,n-k,k-1))]);
  fprintf('w =   %14.5f, p-value = %14.5f\r', [w (1 - cdf('chi2',w,k*(k+1)/2-1))]);
  t = [' Coef' '       Estimate' '           S.E.' ...
      '          Ratio' '     Prob-value'];
  fprintf('%s\r',t);
  d = [linspace(1,k,k)' b se tr pv];
  fprintf('%5d %14.5f %14.5f %14.5f %14.5f\r',d');

%Step 5b: Display graphs.

  plot((1:n),y2,'x')
  xlabel('Index')
  ylabel('e2')
%  scatter(x(:,2),y2,'x')
%  xlabel('X')
%  ylabel('e2')
  
%Step 6: exit to system.

   return;