% Problem Set 14
% Wed. December 5
%Andrew McCallum 2nd edition.

clc; 
clear; 
close all;

%3.
%Input 
%NOTE: All variables are transformed at this stage into logs.
x = log([1:1:15]');
%T = number of observations
T = size(x,1);
%Output
y = log([0.58, 1.10, 1.20, 1.30, 1.95, 2.55, 2.60, 2.90, 3.45, 3.50, 3.60, 4.10, 4.35, 4.40, 4.50]');
%Add constant to Input
x = [ones(15,1), x];
%K = number of regressors
K = size(x,2);
xprimx = x'*x;
xprimy = x'*y;
b = inv(xprimx)*xprimy;
b1 = b(1,1);
b2 = b(2,1);

%3.a
e = y - x*b;
SSE = e'*e;
SSR = (x*b - mean(y,1))'*(x*b - mean(y,1));
SST = (y - mean(y,1))'*(y - mean(y,1));
R2 = SSR/SST; % 0.978597630185609
Fstat = (R2/(K-1))/((1-R2)/(T-K)); % 594.4093716135595
%Testing this using a F(K-1),(T-K) = F(1,13) distribution = 9.07
%We (overwhelmingly) reject the H_0 that all slope coefficients are zero.

%3.b
%H_0: B_2 = 1, H_1: B_2 != 1;
H0 = 1;
sig2 = e'*e/(T-K);

R = [0 1];
r = 1;
%t-test (because we don't know that the innovations are normal we can't
%actually use t-distribution critical values. We must us N(0,1) which for
%two-tailed 10% crit val is 1.64.
tstat = (R*b - r)/(sqrt(sig2)*(R*inv(x'*x)*R')^(0.5)); %-6.618174394971274
% so we (overwhelmingly) reject H0: b2 = 1. 

%Since the restriction is linear we can use the following formula to attain
%the asymtotically robust statistics. Without assuming normality, we cannot
%determine the LR and LM statistics in another way. We'll also use crit vals from
%N(0,1) to test the sig for these. 

F = tstat;

%Wald 
W = F*T/(T-K); %-7.6363550711207

%LR test
LR = T*log(1+F/(T-K)); % -10.672427340715183

%LM test 
LM = (T/(T-K))*(F/(1+F/(T-K))); % -15.55552001395097

%Note/Check: The trinity of tests have the following relationships for their 
%respective %statistics W >= LR >= LM (i.e. the Wald test a priori favors 
%rejection of the %null).

%For one linear restriction we will use a N(0,1) to get the two-tailed critical
%value of 1.64. This applies for each of the three test statistics. 
%Note: If we had J>= 2 restrictions we would use a Chi2(J) distribution.

%All three tests suggest we can overwhelmingly reject the null that b2 = 1.
%I would report the results of each test but recommend using the t-test 
%because it is the simplest, most widely understood and easiest to impliment.

%3.c
%Note: Throughout this question assume normally distributed errors.
r = 0;
tstat = (R*b - r)/(sqrt(sig2)*(R*inv(x'*x)*R')^(0.5)); %24.380512127795978
%Two-tailed t(14) crit val at the 5% level is 2.1448. We overwhelmingly reject 
%H0: b2 = 0.

%Estimate the model under the null H0: b2 = 0;
const = x(:,1);
cprimc = const'*const;
cprimy = const'*y;
b = inv(cprimc)*cprimy;
e = y - const*b;

n = 10000;
mu = 0;
sig2 = e'*e/(T-K);

%Generate values of Y in a world where H0 is true. This will give us b1 and
%sig2 for when b2 = 0.
bmat = b*ones(T,n);
rand('seed',1234);
%Draw new simulated errors from N(0,sig2hat) where sig2hat is the estimate
%given H0 is true.
emat = normrnd(mu,sqrt(sig2),T,n);
yH0mat = bmat + emat;

tstatmat = 9999*ones(1,n);
for i = 1:1:n 
    xprimx = x'*x;
    xprimy = x'*yH0mat(:,i);
    b = inv(xprimx)*xprimy;
    b2 = b(2,1);
    e = yH0mat(:,i) - x*b;
    sig2 = e'*e/(T-K);
    tstatmat(1,i) = (R*b - r)/(sqrt(sig2)*(R*inv(x'*x)*R')^(0.5));
end
     
tstatmat = sort(abs(tstatmat'),'descend');
%Two-tailed 5% critical value for the two-tailed t(14) is 2.1448
pctrejectH0 = max(find(round(tstatmat*100)/100 - 2.1448 > 0))/n; %0.0505

%Generating 10,000 series in a world where H0 is true (giving us b1 and also 
%an estimate of the variance) we find that the calculated
%t-stat is greater than the critical value almost exactly 5% of the time. Hence we
%would reject about 5% of the time. Our empirical size matches the
%theoretical.

%We could determine the empirical finite sample power of the test by
%imposing a given number for b2 taken from a range of values for H1
%(say H1: b2 != 0). We would then do the proceedure
%as we've done above but selecting the empirical t-stats that are less than
%the two-tailed 5% crit val from t(13) (now K = 2 since b2 != 0) which is 2.1604. 
%Doing so would give the frequency of failing to reject given that H0 is false
%(and hence a H1 from our range of values is true) which is just beta. The
% empirical power of the test is then just 1 - beta. The power of the test
% would then be a function over the range of the H1's we've selected to impose.