Contents

AAE 637, Spring 2018, Assignment #5

Charng-Jiun Yu (based on codes from Eduardo Cenci) - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - This code and all accompanying functions are based on previous solutions by: Andrew Schreiber, and Brian Gould (main functions) - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

clear;					    % command to clear memory
clc;                        % command to clear the command window
clear global;               % clear global variables from memory;

Setup

% clear previous log file and start it
delete('question1.txt')
diary('question1.txt');
diary on;

% declare global variables
global numobs do_step critic_limit iter_limit dh func_name depvar rhsvar ...
       numc parnames rhsres parnames_res hetvar hetnames;

% read data
[data,txt,~] = xlsread('RECS_Assign_5_18','sheet1');

data(~any(data(:,3)==1 | data(:,3)==2 | data(:,3)==3, 2),:) = [];
data(~any(data(:,9)==1,2),:) = [];

% change occurence -2 ("not aplicable") in dataset for missing (NaN)
data(data == -2) = NaN;

% extract variables from data
DOEID = data(:,1);
DIVISION = data(:,2);
TYPEHUQ = data(:,3);
HDD30YR = data(:,4);
CDD30YR = data(:,5);
METMIC = data(:,6);
KOWNRENT = data(:,7);
YEARMADE = data(:,8);
HEATHOME = data(:,9);
FUELHEAT = data(:,10);
ADQINSUL = data(:,11);
EDUCATION = data(:,12);
HHAGE = data(:,13);
INCOME = data(:,14);

Generate subsample

NEWENG = (DIVISION==1).*1; % (multiply variables by one to convert logical to double)
MA = (DIVISION==2).*1;
ENC = (DIVISION==3).*1;
WNC = (DIVISION==4).*1;
MOUNT = (DIVISION==8 | DIVISION==9 ).*1;
PACIFIC = (DIVISION==10).*1;
METRO = (METMIC==1).*1;
MICRO = (METMIC==2).*1;
HDD302 = HDD30YR/1000;
CDD302 = CDD30YR/1000;
HDDCDD = HDD302.*CDD302;
HOUSEAGE = 2009-YEARMADE;
OWN = (KOWNRENT==1).*1;
MOBILE = (TYPEHUQ==1).*1;

% define vectors y, X and Z
depvar = (FUELHEAT==3).*1;
rhsvar = horzcat(ones(length(depvar),1),NEWENG,MA,ENC,WNC,PACIFIC,MOUNT,METRO, ...
    MICRO,HDD302,CDD302,HDDCDD,HOUSEAGE,OWN,MOBILE);
hetvar = horzcat(INCOME,HDD302,HOUSEAGE./10,NEWENG);

vars_X = {'CONSTANT','NEWENG', 'MA', 'ENC', 'WNC', 'PACIFIC', 'Mount', 'METRO', 'MICRO',...
          'HDD302', 'CDD302', 'HDDCDD', 'HOUSEAGE', 'OWN', 'MOBILE'};
vars_Z = {'INCOME_het', 'HDD302_het', 'HOUSEAGE_het', 'NEWENG_het'};

% parameter names
parnames = char(vars_X);
hetnames = char(vars_Z);

Question 1(a) - Estimate the probit model

clear vars_X vars_Z;

% Starting values from CRM
b_ini = (rhsvar'*rhsvar)\(rhsvar'*depvar);

% declare funcional form for MLE
func_name = ('probit_llf');

% setttings for MLE:
[numobs,numc] = size(rhsvar);
critic_limit  = 1e-7;
iter_limit    = 500;
do_step       = 1;
dh            = 1e-7;

% type I error probability (significance level for post estimation)
alpha = 0.05;

% run MLE using BHHH for probit model
[b_prob,c_prob,llf] = probit_bwg(b_ini);

% unrestricted LLF
llf_u = llf;
Below is the intermediate iteration information:
  
i =   1LLF = -1963.5720step = 1.0000Paramaters   
 -1.4139  1.0909  0.4056 -0.1325 -0.1177 -0.0412 -0.0857  
 -0.0203 -0.0029  0.0378  0.0207 -0.0092  0.0025  0.0542  
 -0.0550  
  
i =   2LLF = -1420.3023step = 1.0000Paramaters   
 -2.2621  1.4538  0.7023 -0.3140 -0.2983 -0.1228 -0.2335  
 -0.0275  0.0026  0.0869  0.0269 -0.0169  0.0056  0.1170  
 -0.1438  
  
i =   3LLF = -1321.2661step = 1.0000Paramaters   
 -2.7530  1.4534  0.7398 -0.5476 -0.5645 -0.3032 -0.5119  
 -0.0186  0.0252  0.1351 -0.0376 -0.0159  0.0084  0.1897  
 -0.3026  
  
i =   4LLF = -1300.1524step = 1.0000Paramaters   
 -2.7328  1.3147  0.6264 -0.7446 -0.7970 -0.5408 -0.8686  
 -0.0155  0.0388  0.1433 -0.2369  0.0072  0.0096  0.2337  
 -0.4682  
  
i =   5LLF = -1294.0470step = 1.0000Paramaters   
 -2.5751  1.2191  0.5352 -0.8609 -0.9320 -0.6276 -1.0785  
 -0.0379  0.0263  0.1289 -0.4956  0.0594  0.0097  0.2460  
 -0.5330  
  
i =   6LLF = -1292.0712step = 1.0000Paramaters   
 -2.4945  1.1672  0.4809 -0.9370 -1.0281 -0.6123 -1.1440  
 -0.0647  0.0114  0.1171 -0.7344  0.1189  0.0097  0.2494  
 -0.5437  
  
i =   7LLF = -1291.8383step = 1.0000Paramaters   
 -2.4619  1.1466  0.4599 -0.9688 -1.0692 -0.6037 -1.1640  
 -0.0757  0.0060  0.1118 -0.8479  0.1464  0.0096  0.2506  
 -0.5463  
  
i =   8LLF = -1291.8354step = 1.0000Paramaters   
 -2.4577  1.1442  0.4576 -0.9724 -1.0739 -0.6031 -1.1663  
 -0.0768  0.0055  0.1112 -0.8620  0.1497  0.0096  0.2508  
 -0.5466  
  
i =   9LLF = -1291.8354step = 1.0000Paramaters   
 -2.4577  1.1442  0.4576 -0.9724 -1.0739 -0.6031 -1.1663  
 -0.0769  0.0055  0.1112 -0.8622  0.1497  0.0096  0.2508  
 -0.5466  
  
i =  10LLF = -1291.8354step = 0.5000Paramaters   
 -2.4577  1.1442  0.4576 -0.9724 -1.0739 -0.6031 -1.1663  
 -0.0769  0.0055  0.1112 -0.8622  0.1497  0.0096  0.2508  
 -0.5466  
  
i =  11LLF = -1291.8354step = 1.0000Paramaters   
 -2.4577  1.1442  0.4576 -0.9724 -1.0739 -0.6031 -1.1663  
 -0.0769  0.0055  0.1112 -0.8622  0.1497  0.0096  0.2508  
 -0.5466  
   
This is the Probit Maximum Likelihood Result
---------------------------------------------------------- 
 Variables      Value    Std.Err    Z-Value    P-Value  
---------------------------------------------------------- 
  CONSTANT   -2.45767    0.41898   -5.86585    0.00000  
  NEWENG      1.14421    0.12034    9.50846    0.00000  
  MA          0.45760    0.10941    4.18245    0.00003  
  ENC        -0.97243    0.14845   -6.55042    0.00000  
  WNC        -1.07392    0.16624   -6.46020    0.00000  
  PACIFIC    -0.60314    0.20898   -2.88606    0.00390  
  Mount      -1.16634    0.23438   -4.97632    0.00000  
  METRO      -0.07686    0.11460   -0.67064    0.50245  
  MICRO       0.00549    0.13836    0.03968    0.96835  
  HDD302      0.11118    0.04451    2.49792    0.01249  
  CDD302     -0.86218    0.21663   -3.97990    0.00007  
  HDDCDD      0.14975    0.05444    2.75088    0.00594  
  HOUSEAGE    0.00962    0.00107    9.00750    0.00000  
  OWN         0.25078    0.08542    2.93588    0.00333  
  MOBILE     -0.54662    0.19543   -2.79694    0.00516  
  

Test if exogenous variables are jointly significant by estimating the restricted model

% (1) First way: estimate the restricted model

numc_full = numc; % store number of parameters for the full model

func_name = ('probit_llf_res');
rhsres = rhsvar(:,1);
[~,nl] = size(parnames);
parnames_res = parnames(1, 1:nl);
[~,numc] = size(rhsres);

% initial values for betas
b_ini_res = (rhsres'*rhsres)\(rhsres'*depvar);

% run MLE using BHHH for restricted probit model
[~,~,llf] = probit_bwg_res(b_ini_res);

% restricted LLF
llf_r = sum(llf);

% perform likelihood ratio test
l_ratio = 2*(llf_u - llf_r);
fprintf('\nThe likelihood ratio value is: %4.4f\n', l_ratio);
crit = chi_bwg(alpha,numc_full-numc,l_ratio);

% (2) Second way: Use the formula from lecture slide "qualitative_2_12",
% page 28

ybar = mean(depvar);
llf_r = numobs*(ybar*log(ybar)+(1-ybar)*log(1-ybar));

% perform likelihood ratio test
l_ratio = 2*(llf_u - llf_r);
fprintf('\nThe likelihood ratio value is: %4.4f\n', l_ratio);
crit = chi_bwg(alpha,numc_full-numc,l_ratio);
Below is the intermediate iteration information:
  
i =   1LLF = -2502.7118step = 1.0000Paramaters  -1.0598  
  
  
i =   2LLF = -2253.9277step = 1.0000Paramaters  -1.4243  
  
  
i =   3LLF = -2249.0528step = 1.0000Paramaters  -1.4850  
  
  
i =   4LLF = -2249.0494step = 1.0000Paramaters  -1.4867  
  
  
i =   5LLF = -2249.0494step = 1.0000Paramaters  -1.4867  
  
  
i =   6LLF = -2249.0494step = 1.0000Paramaters  -1.4867  
  
   
This is the Probit Maximum Likelihood Result
---------------------------------------------------------- 
 Variables      Value    Std.Err    Z-Value    P-Value  
---------------------------------------------------------- 
  CONSTANT   -1.48669    0.02016  -73.74766    0.00000  
  

The likelihood ratio value is: 1914.4281
This is the Prob. of Type I error:  0.05  
This is the critical Chi-Sq. Stat:  23.685-->Reject H0 
The likelihood ratio value is: 1914.4281
This is the Prob. of Type I error:  0.05  
This is the critical Chi-Sq. Stat:  23.685-->Reject H0 

Question 1(b).i - Test the impact of geographic locations:

NewEng, MA, ENC, WNC, Pacific, and Mount (columns 2-7)

func_name = ('probit_llf_res');
rhsres = rhsvar(:,[1 8:15]);
[~,nl] = size(parnames);
parnames_res = parnames([1 8:15], 1:nl);
[numobs,numc] = size(rhsres);

% initial values for betas
b_ini_res = (rhsres'*rhsres)\(rhsres'*depvar);

% run MLE using BHHH for restricted probit model
[~,~,llf] = probit_bwg_res(b_ini_res);

% restricted LLF
llf_r = sum(llf);

% perform likelihood ratio test
l_ratio = 2*(llf_u - llf_r);
fprintf('\nThe likelihood ratio value is: %4.4f\n', l_ratio);
crit = chi_bwg(alpha,numc_full-numc,l_ratio);

% re-initialize number of obs and parameters
[numobs,numc] = size(rhsvar);
Below is the intermediate iteration information:
  
i =   1LLF = -2302.8564step = 1.0000Paramaters   
 -1.4968  0.0396 -0.0097  0.0842  0.0351 -0.0489  0.0039  
  0.0714 -0.0680  
  
i =   2LLF = -1874.2376step = 1.0000Paramaters   
 -2.3783  0.1302  0.0152  0.1602  0.0236 -0.0902  0.0080  
  0.1605 -0.1749  
  
i =   3LLF = -1806.1045step = 1.0000Paramaters   
 -2.7468  0.2228  0.0518  0.1855 -0.1018 -0.0953  0.0106  
  0.2356 -0.3388  
  
i =   4LLF = -1792.0462step = 1.0000Paramaters   
 -2.6952  0.2377  0.0501  0.1695 -0.3068 -0.0649  0.0112  
  0.2615 -0.4528  
  
i =   5LLF = -1787.9630step = 1.0000Paramaters   
 -2.5424  0.2203  0.0336  0.1456 -0.5173 -0.0255  0.0112  
  0.2652 -0.4709  
  
i =   6LLF = -1787.3221step = 1.0000Paramaters   
 -2.4572  0.2103  0.0256  0.1324 -0.6397 -0.0025  0.0112  
  0.2660 -0.4685  
  
i =   7LLF = -1787.3071step = 1.0000Paramaters   
 -2.4435  0.2088  0.0245  0.1302 -0.6615  0.0015  0.0112  
  0.2662 -0.4680  
  
i =   8LLF = -1787.3071step = 1.0000Paramaters   
 -2.4433  0.2088  0.0245  0.1302 -0.6620  0.0016  0.0112  
  0.2662 -0.4679  
  
i =   9LLF = -1787.3071step = 1.0000Paramaters   
 -2.4433  0.2088  0.0245  0.1302 -0.6620  0.0016  0.0112  
  0.2662 -0.4679  
  
i =  10LLF = -1787.3071step = 0.2500Paramaters   
 -2.4433  0.2088  0.0245  0.1302 -0.6620  0.0016  0.0112  
  0.2662 -0.4679  
   
This is the Probit Maximum Likelihood Result
---------------------------------------------------------- 
 Variables      Value    Std.Err    Z-Value    P-Value  
---------------------------------------------------------- 
  CONSTANT   -2.44325    0.20824  -11.73258    0.00000  
  METRO       0.20880    0.09528    2.19153    0.02841  
  MICRO       0.02451    0.11747    0.20866    0.83472  
  HDD302      0.13019    0.02329    5.59050    0.00000  
  CDD302     -0.66204    0.13841   -4.78322    0.00000  
  HDDCDD      0.00161    0.02876    0.05585    0.95546  
  HOUSEAGE    0.01120    0.00089   12.52381    0.00000  
  OWN         0.26616    0.07299    3.64681    0.00027  
  MOBILE     -0.46794    0.15951   -2.93356    0.00335  
  

The likelihood ratio value is: 990.9434
This is the Prob. of Type I error:  0.05  
This is the critical Chi-Sq. Stat:  12.592-->Reject H0 

Question 1(b).ii - Calculate the discrete change effect for Mobile

diff = disc_mobile(b_prob);
fprintf('\nThe average discrete change effect for dummy Mobile is: %8.4f\n', diff);

% Delta Method to obtain variance of discrete change effect
diff_der = Grad(b_prob,'disc_mobile',1);
diff_cov = diff_der * c_prob * diff_der';

% compute the Z-value
zstat = (diff - 0)/sqrt(diff_cov);
fprintf('\nThe z-stat for testing discrete change = 0 is %8.4f \n',...
        zstat);
The average discrete change effect for dummy Mobile is:  -0.0342

The z-stat for testing discrete change = 0 is  -3.5945 

Question 1(b).iii - Calculate and test the marginal effect of HDD302

diff = marg_hdd(b_prob);
fprintf('\nThe average marginal effect of HDD302 is: %8.4f\n', diff);

% Delta Method to obtain variance of discrete change effect
diff_der = Grad(b_prob,'marg_hdd',1);
diff_cov = diff_der * c_prob * diff_der';

% perform t-test
tstat = (diff - 0)/sqrt(diff_cov);
pvalue = 2*(1-tcdf(abs(tstat),numobs-numc));
fprintf('\nThe t-stat for testing marginal effect of HDD302 = 0 is %8.4f (p-value = %5.4f)\n',...
        tstat, pvalue);
if pvalue<0.05
    disp('Reject H_0 under 0.05 significance level');
else
    disp('Fail to reject H_0 under 0.05 significance level');
end
The average marginal effect of HDD302 is:   0.0172

The t-stat for testing marginal effect of HDD302 = 0 is   3.6792 (p-value = 0.0002)
Reject H_0 under 0.05 significance level

Question 1(b).iv - Test the significance of average interaction effect of HDD302 and CDD302

diff = iter_effect(b_prob);
fprintf('\nThe average interaction effect of HDD302 and CDD302 is: %8.4f\n', diff);

% Delta Method to obtain variance of discrete change effect
diff_der = Grad(b_prob,'iter_effect',1);
diff_cov = diff_der * c_prob * diff_der';

% perform t-test
tstat = (diff - 0)/sqrt(diff_cov);
pvalue = 2*(1-tcdf(abs(tstat),numobs-numc));
fprintf('\nThe t-stat for testing interaction effect of HDD302 and CDD302 = 0 is %8.4f (p-value = %5.4f)\n',...
        tstat, pvalue);
if pvalue<0.05
    disp('Reject H_0 under 0.05 significance level');
else
    disp('Fail to reject H_0 under 0.05 significance level');
end
The average interaction effect of HDD302 and CDD302 is:   0.0101

The t-stat for testing interaction effect of HDD302 and CDD302 = 0 is   1.5401 (p-value = 0.1236)
Fail to reject H_0 under 0.05 significance level

Question 1(c) - Calculate the average elasticity of HDD302

diff = elast_hdd(b_prob);
fprintf('\nThe average elasticity of HDD302 is: %8.4f\n', diff);

% Delta Method to obtain variance of discrete change effect
diff_der = Grad(b_prob,'elast_hdd',1);
diff_cov = diff_der * c_prob * diff_der';

% perform t-test
tstat = (diff - 0)/sqrt(diff_cov);
pvalue = 2*(1-tcdf(abs(tstat),numobs-numc));
fprintf('\nThe t-stat for testing average elasticity of HDD302 = 0 is %8.4f (p-value = %5.4f)\n',...
        tstat, pvalue);
if pvalue<0.05
    disp('Reject H_0 under 0.05 significance level');
else
    disp('Fail to reject H_0 under 0.05 significance level');
end
The average elasticity of HDD302 is:   2.8774

The t-stat for testing average elasticity of HDD302 = 0 is   3.5322 (p-value = 0.0004)
Reject H_0 under 0.05 significance level

Question 1(d) - Calculate the average elasticity of HouseAge

e_avg = elast_hage_avg(b_prob);
fprintf('\nThe average elasticity of HouseAge is: %8.4f\n', e_avg);

% calculate the elasticity of HouseAge at the mean of the data
e_mean = elast_hage_mean(b_prob);
fprintf('\nThe elasticity of HouseAge at the mean of the data is: %8.4f\n', e_mean);

% calculate their difference
diff = elast_hage_diff(b_prob);
fprintf('\nThe difference of these two elasticities is: %8.4f\n', diff);

% Delta Method to obtain variance of discrete change effect
diff_der = Grad(b_prob,'elast_hage_diff',1);
diff_cov = diff_der * c_prob * diff_der';

% perform t-test
tstat = (diff - 0)/sqrt(diff_cov);
pvalue = 2*(1-tcdf(abs(tstat),numobs-numc));
fprintf('\nThe t-stat for testing equal elasticities is %8.4f (p-value = %5.4f)\n',...
        tstat, pvalue);
if pvalue<0.05
    disp('Reject H_0 under 0.05 significance level');
else
    disp('Fail to reject H_0 under 0.05 significance level');
end

clear b_ini a_ini p_ini;
diary off;
The average elasticity of HouseAge is:   0.8900

The elasticity of HouseAge at the mean of the data is:   0.9768

The difference of these two elasticities is:  -0.0868

The t-stat for testing equal elasticities is  -5.7188 (p-value = 0.0000)
Reject H_0 under 0.05 significance level

Question 2(a) - Estimate heteroskedastic model

clear previous log file and start it

delete('question2.txt')
diary('question2.txt');
diary on;

% starting values
b_ini = b_prob;
g_ini = [.01, .01, .01, 0.01]';
p_ini = vertcat(b_ini, g_ini);

names = char(parnames,hetnames);

% declare funcional form for MLE
func_name = ('probit_hetero_llf');

% setttings for MLE:
[numobs,numc] = size(rhsvar);
critic_limit  = 1e-6;
iter_limit    = 500;
do_step       = 1;
dh            = 1e-7;

[b_phet,c_phet,llf_het] = max_bhhh(p_ini,names);

% unrestricted LLF (het.)
llf_het_u = sum(llf_het);

% perform likelihood ratio test
l_ratio = 2*(llf_het_u - llf_u);
fprintf('\nThe likelihood ratio value is: %4.4f\n', l_ratio);
crit = chi_bwg(alpha,length(b_phet)-length(b_prob),l_ratio);
disp('');
Below is the intermediate iteration information:
  
  
  
i =   1LLF = -1297.5748step = 0.250Parameters  -3.2485  1.5419  0.6267 -1.3381 -1.5406 -0.6471 -1.4273 -0.1377 -0.0384  0.1571 -0.7475  0.1409  0.0112  0.3353  
 -0.7369 -0.0335  0.0497  0.0241 -0.0664  
  
i =   2LLF = -1287.1356step = 1.000Parameters  -3.8397  2.2451  0.9425 -1.9436 -2.2465 -0.6805 -2.0034 -0.2178 -0.1359  0.1474 -0.7875  0.1343  0.0111  0.4923  
 -1.0689 -0.0682  0.0779  0.0488 -0.2070  
  
i =   3LLF = -1285.5030step = 0.500Parameters  -3.3447  1.9413  0.7933 -1.6680 -1.8562 -0.7224 -1.8913 -0.1038  0.0067  0.1117 -1.0072  0.1716  0.0094  0.4615  
 -0.8194 -0.0493  0.0605  0.0402 -0.2326  
  
i =   4LLF = -1284.8400step = 0.500Parameters  -3.5793  2.1387  0.8846 -1.7845 -2.0034 -0.7860 -1.9632 -0.1214 -0.0128  0.1213 -0.8467  0.1368  0.0093  0.4944  
 -0.8002 -0.0457  0.0684  0.0450 -0.3295  
  
i =   5LLF = -1284.3227step = 1.000Parameters  -3.3337  2.0526  0.8492 -1.6987 -1.8629 -0.7976 -1.9001 -0.0714  0.0573  0.1036 -0.9458  0.1448  0.0079  0.4736  
 -0.6442 -0.0348  0.0585  0.0461 -0.4489  
  
i =   6LLF = -1284.1521step = 0.500Parameters  -3.5449  2.2118  0.9214 -1.8121 -2.0113 -0.8462 -1.9762 -0.0776  0.0569  0.1123 -0.8429  0.1241  0.0079  0.4946  
 -0.6407 -0.0367  0.0673  0.0497 -0.5378  
  
i =   7LLF = -1283.9810step = 1.000Parameters  -3.3937  2.1866  0.9146 -1.7773 -1.9460 -0.8593 -1.9697 -0.0444  0.1017  0.0991 -0.9132  0.1275  0.0071  0.4744  
 -0.5674 -0.0320  0.0637  0.0514 -0.6396  
  
i =   8LLF = -1283.9019step = 0.500Parameters  -3.5233  2.2915  0.9610 -1.8558 -2.0493 -0.8996 -2.0269 -0.0449  0.1063  0.1037 -0.8464  0.1134  0.0071  0.4845  
 -0.5589 -0.0331  0.0697  0.0532 -0.6983  
  
i =   9LLF = -1283.8182step = 1.000Parameters  -3.4193  2.2877  0.9628 -1.8350 -2.0070 -0.9121 -2.0298 -0.0237  0.1341  0.0933 -0.8927  0.1149  0.0066  0.4676  
 -0.5100 -0.0297  0.0680  0.0545 -0.7714  
  
i =  10LLF = -1283.7803step = 0.500Parameters  -3.5036  2.3591  0.9936 -1.8911 -2.0817 -0.9416 -2.0715 -0.0228  0.1386  0.0959 -0.8494  0.1057  0.0065  0.4723  
 -0.5023 -0.0305  0.0721  0.0555 -0.8115  
  
i =  11LLF = -1283.7394step = 1.000Parameters  -3.4350  2.3624  0.9976 -1.8789 -2.0550 -0.9508 -2.0761 -0.0091  0.1559  0.0885 -0.8803  0.1066  0.0062  0.4595  
 -0.4697 -0.0280  0.0712  0.0565 -0.8628  
  
i =  12LLF = -1283.7213step = 0.500Parameters  -3.4897  2.4109  1.0182 -1.9182 -2.1075 -0.9718 -2.1058 -0.0080  0.1594  0.0900 -0.8521  0.1005  0.0062  0.4618  
 -0.4636 -0.0286  0.0740  0.0571 -0.8903  
  
i =  13LLF = -1283.7017step = 1.000Parameters  -3.4451  2.4159  1.0221 -1.9112 -2.0910 -0.9785 -2.1104  0.0008  0.1703  0.0848 -0.8724  0.1011  0.0059  0.4527  
 -0.4419 -0.0268  0.0735  0.0578 -0.9256  
  
i =  14LLF = -1283.6934step = 0.500Parameters  -3.4803  2.4485  1.0358 -1.9380 -2.1268 -0.9931 -2.1309  0.0017  0.1727  0.0856 -0.8540  0.0971  0.0059  0.4537  
 -0.4374 -0.0272  0.0754  0.0581 -0.9443  
  
i =  15LLF = -1283.6842step = 1.000Parameters  -3.4518  2.4535  1.0391 -1.9343 -2.1170 -0.9979 -2.1349  0.0073  0.1796  0.0821 -0.8673  0.0974  0.0058  0.4474  
 -0.4229 -0.0260  0.0751  0.0586 -0.9684  
  
i =  16LLF = -1283.6802step = 1.000Parameters  -3.4966  2.4966  1.0572 -1.9698 -2.1646 -1.0176 -2.1624  0.0087  0.1830  0.0831 -0.8436  0.0922  0.0058  0.4483  
 -0.4164 -0.0264  0.0776  0.0591 -0.9935  
  
i =  17LLF = -1283.6762step = 0.500Parameters  -3.4564  2.4793  1.0507 -1.9505 -2.1353 -1.0112 -2.1521  0.0116  0.1857  0.0803 -0.8641  0.0950  0.0057  0.4437  
 -0.4102 -0.0254  0.0763  0.0591 -0.9972  
  
i =  18LLF = -1283.6743step = 1.000Parameters  -3.4845  2.5074  1.0625 -1.9733 -2.1659 -1.0242 -2.1700  0.0127  0.1880  0.0808 -0.8486  0.0915  0.0057  0.4440  
 -0.4057 -0.0256  0.0779  0.0594 -1.0137  
  
i =  19LLF = -1283.6726step = 0.500Parameters  -3.4588  2.4965  1.0584 -1.9612 -2.1475 -1.0202 -2.1636  0.0145  0.1897  0.0790 -0.8620  0.0934  0.0056  0.4411  
 -0.4017 -0.0250  0.0770  0.0595 -1.0164  
  
i =  20LLF = -1283.6717step = 1.000Parameters  -3.4768  2.5148  1.0661 -1.9759 -2.1671 -1.0287 -2.1752  0.0153  0.1912  0.0793 -0.8520  0.0911  0.0056  0.4411  
 -0.3985 -0.0251  0.0780  0.0597 -1.0273  
  
i =  21LLF = -1283.6710step = 0.500Parameters  -3.4605  2.5081  1.0636 -1.9684 -2.1557 -1.0262 -2.1713  0.0164  0.1923  0.0781 -0.8606  0.0924  0.0056  0.4392  
 -0.3959 -0.0247  0.0775  0.0597 -1.0292  
  
i =  22LLF = -1283.6706step = 1.000Parameters  -3.4717  2.5199  1.0686 -1.9778 -2.1681 -1.0318 -2.1788  0.0170  0.1933  0.0783 -0.8543  0.0909  0.0056  0.4392  
 -0.3938 -0.0248  0.0782  0.0599 -1.0364  
  
i =  23LLF = -1283.6703step = 0.500Parameters  -3.4615  2.5158  1.0670 -1.9732 -2.1612 -1.0303 -2.1765  0.0177  0.1940  0.0775 -0.8598  0.0917  0.0055  0.4380  
 -0.3921 -0.0245  0.0779  0.0599 -1.0377  
  
i =  24LLF = -1283.6701step = 1.000Parameters  -3.4685  2.5234  1.0703 -1.9792 -2.1690 -1.0339 -2.1813  0.0181  0.1947  0.0776 -0.8557  0.0907  0.0055  0.4379  
 -0.3907 -0.0245  0.0783  0.0600 -1.0424  
  
i =  25LLF = -1283.6700step = 0.500Parameters  -3.4621  2.5209  1.0693 -1.9764 -2.1648 -1.0330 -2.1800  0.0185  0.1952  0.0771 -0.8592  0.0912  0.0055  0.4372  
 -0.3896 -0.0244  0.0781  0.0600 -1.0434  
  
i =  26LLF = -1283.6699step = 1.000Parameters  -3.4665  2.5258  1.0714 -1.9802 -2.1697 -1.0353 -2.1830  0.0188  0.1956  0.0772 -0.8567  0.0906  0.0055  0.4371  
 -0.3886 -0.0244  0.0783  0.0601 -1.0464  
  
i =  27LLF = -1283.6698step = 0.500Parameters  -3.4625  2.5243  1.0708 -1.9785 -2.1672 -1.0348 -2.1823  0.0191  0.1959  0.0769 -0.8589  0.0909  0.0055  0.4366  
 -0.3879 -0.0243  0.0782  0.0601 -1.0471  
  
i =  28LLF = -1283.6698step = 1.000Parameters  -3.4652  2.5274  1.0722 -1.9809 -2.1702 -1.0363 -2.1842  0.0193  0.1962  0.0769 -0.8573  0.0905  0.0055  0.4366  
 -0.3872 -0.0243  0.0784  0.0601 -1.0491  
  
i =  29LLF = -1283.6698step = 0.500Parameters  -3.4627  2.5265  1.0718 -1.9799 -2.1687 -1.0360 -2.1838  0.0194  0.1964  0.0767 -0.8587  0.0907  0.0055  0.4363  
 -0.3868 -0.0242  0.0783  0.0601 -1.0496  
  
i =  30LLF = -1283.6697step = 1.000Parameters  -3.4644  2.5285  1.0727 -1.9814 -2.1706 -1.0369 -2.1850  0.0196  0.1966  0.0767 -0.8577  0.0904  0.0055  0.4362  
 -0.3863 -0.0242  0.0784  0.0601 -1.0509  
  
i =  31LLF = -1283.6697step = 0.500Parameters  -3.4629  2.5280  1.0725 -1.9808 -2.1697 -1.0368 -2.1848  0.0197  0.1967  0.0766 -0.8585  0.0906  0.0055  0.4360  
 -0.3861 -0.0242  0.0784  0.0601 -1.0512  
  
i =  32LLF = -1283.6697step = 1.000Parameters  -3.4639  2.5293  1.0730 -1.9817 -2.1709 -1.0374 -2.1856  0.0198  0.1969  0.0766 -0.8579  0.0904  0.0055  0.4360  
 -0.3858 -0.0242  0.0784  0.0602 -1.0521  
  
i =  33LLF = -1283.6697step = 0.500Parameters  -3.4629  2.5289  1.0729 -1.9814 -2.1704 -1.0373 -2.1854  0.0198  0.1969  0.0765 -0.8584  0.0905  0.0055  0.4359  
 -0.3856 -0.0242  0.0784  0.0602 -1.0523  
  
i =  34LLF = -1283.6697step = 1.000Parameters  -3.4636  2.5298  1.0733 -1.9820 -2.1711 -1.0377 -2.1859  0.0199  0.1970  0.0765 -0.8580  0.0904  0.0055  0.4358  
 -0.3854 -0.0242  0.0785  0.0602 -1.0528  
  
i =  35LLF = -1283.6697step = 0.500Parameters  -3.4630  2.5296  1.0732 -1.9818 -2.1708 -1.0376 -2.1859  0.0199  0.1971  0.0764 -0.8584  0.0904  0.0055  0.4357  
 -0.3852 -0.0242  0.0784  0.0602 -1.0530  
  
i =  36LLF = -1283.6697step = 1.000Parameters  -3.4634  2.5301  1.0734 -1.9821 -2.1713 -1.0379 -2.1862  0.0200  0.1971  0.0764 -0.8581  0.0904  0.0055  0.4357  
 -0.3851 -0.0242  0.0785  0.0602 -1.0533  
  
i =  37LLF = -1283.6697step = 0.500Parameters  -3.4630  2.5300  1.0734 -1.9820 -2.1711 -1.0378 -2.1861  0.0200  0.1972  0.0764 -0.8584  0.0904  0.0055  0.4357  
 -0.3850 -0.0241  0.0785  0.0602 -1.0534  
  
i =  38LLF = -1283.6697step = 1.000Parameters  -3.4633  2.5303  1.0735 -1.9822 -2.1714 -1.0380 -2.1863  0.0200  0.1972  0.0764 -0.8582  0.0903  0.0055  0.4357  
 -0.3849 -0.0241  0.0785  0.0602 -1.0537  
  
i =  39LLF = -1283.6697step = 0.500Parameters  -3.4630  2.5303  1.0735 -1.9822 -2.1713 -1.0380 -2.1863  0.0200  0.1972  0.0764 -0.8583  0.0904  0.0055  0.4356  
 -0.3849 -0.0241  0.0785  0.0602 -1.0538  
  
i =  40LLF = -1283.6697step = 1.000Parameters  -3.4632  2.5305  1.0736 -1.9823 -2.1714 -1.0381 -2.1865  0.0201  0.1973  0.0764 -0.8582  0.0903  0.0055  0.4356  
 -0.3848 -0.0241  0.0785  0.0602 -1.0539  
  
i =  41LLF = -1283.6697step = 0.500Parameters  -3.4630  2.5304  1.0736 -1.9823 -2.1714 -1.0381 -2.1864  0.0201  0.1973  0.0764 -0.8583  0.0904  0.0055  0.4356  
 -0.3848 -0.0241  0.0785  0.0602 -1.0540  
  
i =  42LLF = -1283.6697step = 1.000Parameters  -3.4631  2.5306  1.0736 -1.9824 -2.1715 -1.0381 -2.1865  0.0201  0.1973  0.0764 -0.8583  0.0903  0.0055  0.4356  
 -0.3848 -0.0241  0.0785  0.0602 -1.0541  
  
i =  43LLF = -1283.6697step = 1.000Parameters  -3.4629  2.5305  1.0736 -1.9823 -2.1714 -1.0381 -2.1865  0.0201  0.1973  0.0763 -0.8584  0.0903  0.0055  0.4356  
 -0.3847 -0.0241  0.0785  0.0602 -1.0541  
  
i =  44LLF = -1283.6697step = 0.500Parameters  -3.4631  2.5307  1.0737 -1.9824 -2.1715 -1.0382 -2.1866  0.0201  0.1973  0.0763 -0.8583  0.0903  0.0055  0.4356  
 -0.3847 -0.0241  0.0785  0.0602 -1.0542  
  
i =  45LLF = -1283.6697step = 1.000Parameters  -3.4630  2.5306  1.0736 -1.9824 -2.1715 -1.0382 -2.1866  0.0201  0.1973  0.0763 -0.8584  0.0903  0.0055  0.4356  
 -0.3847 -0.0241  0.0785  0.0602 -1.0542  
  
i =  46LLF = -1283.6697step = 0.500Parameters  -3.4631  2.5307  1.0737 -1.9824 -2.1715 -1.0382 -2.1866  0.0201  0.1973  0.0763 -0.8583  0.0903  0.0055  0.4356  
 -0.3847 -0.0241  0.0785  0.0602 -1.0542  
  
i =  47LLF = -1283.6697step = 1.000Parameters  -3.4630  2.5307  1.0737 -1.9824 -2.1715 -1.0382 -2.1866  0.0201  0.1973  0.0763 -0.8583  0.0903  0.0055  0.4355  
 -0.3847 -0.0241  0.0785  0.0602 -1.0542  
  
i =  48LLF = -1283.6697step = 0.500Parameters  -3.4631  2.5307  1.0737 -1.9825 -2.1716 -1.0382 -2.1866  0.0201  0.1973  0.0763 -0.8583  0.0903  0.0055  0.4355  
 -0.3847 -0.0241  0.0785  0.0602 -1.0543  
  
i =  49LLF = -1283.6697step = 1.000Parameters  -3.4630  2.5307  1.0737 -1.9824 -2.1715 -1.0382 -2.1866  0.0201  0.1973  0.0763 -0.8583  0.0903  0.0055  0.4355  
 -0.3846 -0.0241  0.0785  0.0602 -1.0543  
  
i =  50LLF = -1283.6697step = 0.500Parameters  -3.4631  2.5307  1.0737 -1.9825 -2.1716 -1.0382 -2.1866  0.0201  0.1973  0.0763 -0.8583  0.0903  0.0055  0.4355  
 -0.3846 -0.0241  0.0785  0.0602 -1.0543  
  
i =  51LLF = -1283.6697step = 1.000Parameters  -3.4630  2.5307  1.0737 -1.9825 -2.1716 -1.0382 -2.1866  0.0201  0.1973  0.0763 -0.8583  0.0903  0.0055  0.4355  
 -0.3846 -0.0241  0.0785  0.0602 -1.0543  
  
i =  52LLF = -1283.6697step = 0.500Parameters  -3.4630  2.5308  1.0737 -1.9825 -2.1716 -1.0382 -2.1867  0.0201  0.1973  0.0763 -0.8583  0.0903  0.0055  0.4355  
 -0.3846 -0.0241  0.0785  0.0602 -1.0543  
  
i =  53LLF = -1283.6697step = 1.000Parameters  -3.4630  2.5308  1.0737 -1.9825 -2.1716 -1.0383 -2.1867  0.0201  0.1973  0.0763 -0.8583  0.0903  0.0055  0.4355  
 -0.3846 -0.0241  0.0785  0.0602 -1.0543  
  
i =  54LLF = -1283.6697step = 0.500Parameters  -3.4630  2.5308  1.0737 -1.9825 -2.1716 -1.0383 -2.1867  0.0201  0.1974  0.0763 -0.8583  0.0903  0.0055  0.4355  
 -0.3846 -0.0241  0.0785  0.0602 -1.0543  
  
i =  55LLF = -1283.6697step = 1.000Parameters  -3.4630  2.5308  1.0737 -1.9825 -2.1716 -1.0383 -2.1867  0.0201  0.1974  0.0763 -0.8583  0.0903  0.0055  0.4355  
 -0.3846 -0.0241  0.0785  0.0602 -1.0543  
  
i =  56LLF = -1283.6697step = 0.500Parameters  -3.4630  2.5308  1.0737 -1.9825 -2.1716 -1.0383 -2.1867  0.0201  0.1974  0.0763 -0.8583  0.0903  0.0055  0.4355  
 -0.3846 -0.0241  0.0785  0.0602 -1.0543  
  
i =  57LLF = -1283.6697step = 1.000Parameters  -3.4630  2.5308  1.0737 -1.9825 -2.1716 -1.0383 -2.1867  0.0201  0.1974  0.0763 -0.8583  0.0903  0.0055  0.4355  
 -0.3846 -0.0241  0.0785  0.0602 -1.0543  
  
i =  58LLF = -1283.6697step = 0.500Parameters  -3.4630  2.5308  1.0737 -1.9825 -2.1716 -1.0383 -2.1867  0.0201  0.1974  0.0763 -0.8583  0.0903  0.0055  0.4355  
 -0.3846 -0.0241  0.0785  0.0602 -1.0543  
  
i =  59LLF = -1283.6697step = 1.000Parameters  -3.4630  2.5308  1.0737 -1.9825 -2.1716 -1.0383 -2.1867  0.0201  0.1974  0.0763 -0.8583  0.0903  0.0055  0.4355  
 -0.3846 -0.0241  0.0785  0.0602 -1.0543  
  
i =  60LLF = -1283.6697step = 0.500Parameters  -3.4630  2.5308  1.0737 -1.9825 -2.1716 -1.0383 -2.1867  0.0201  0.1974  0.0763 -0.8583  0.0903  0.0055  0.4355  
 -0.3846 -0.0241  0.0785  0.0602 -1.0543  
  
i =  61LLF = -1283.6697step = 1.000Parameters  -3.4630  2.5308  1.0737 -1.9825 -2.1716 -1.0383 -2.1867  0.0201  0.1974  0.0763 -0.8583  0.0903  0.0055  0.4355  
 -0.3846 -0.0241  0.0785  0.0602 -1.0543  
  
i =  62LLF = -1283.6697step = 0.500Parameters  -3.4630  2.5308  1.0737 -1.9825 -2.1716 -1.0383 -2.1867  0.0201  0.1974  0.0763 -0.8583  0.0903  0.0055  0.4355  
 -0.3846 -0.0241  0.0785  0.0602 -1.0543  
  
i =  63LLF = -1283.6697step = 0.500Parameters  -3.4630  2.5308  1.0737 -1.9825 -2.1716 -1.0383 -2.1867  0.0201  0.1974  0.0763 -0.8583  0.0903  0.0055  0.4355  
 -0.3846 -0.0241  0.0785  0.0602 -1.0543  
  
i =  64LLF = -1283.6697step = 1.000Parameters  -3.4630  2.5308  1.0737 -1.9825 -2.1716 -1.0383 -2.1867  0.0201  0.1974  0.0763 -0.8583  0.0903  0.0055  0.4355  
 -0.3846 -0.0241  0.0785  0.0602 -1.0543  
  
i =  65LLF = -1283.6697step = 0.500Parameters  -3.4630  2.5308  1.0737 -1.9825 -2.1716 -1.0383 -2.1867  0.0201  0.1974  0.0763 -0.8583  0.0903  0.0055  0.4355  
 -0.3846 -0.0241  0.0785  0.0602 -1.0543  
  
i =  66LLF = -1283.6697step = 1.000Parameters  -3.4630  2.5308  1.0737 -1.9825 -2.1716 -1.0383 -2.1867  0.0201  0.1974  0.0763 -0.8583  0.0903  0.0055  0.4355  
 -0.3846 -0.0241  0.0785  0.0602 -1.0543  
  
i =  67LLF = -1283.6697step = 0.500Parameters  -3.4630  2.5308  1.0737 -1.9825 -2.1716 -1.0383 -2.1867  0.0201  0.1974  0.0763 -0.8583  0.0903  0.0055  0.4355  
 -0.3846 -0.0241  0.0785  0.0602 -1.0543  
  
i =  68LLF = -1283.6697step = 0.250Parameters  -3.4630  2.5308  1.0737 -1.9825 -2.1716 -1.0383 -2.1867  0.0201  0.1974  0.0763 -0.8583  0.0903  0.0055  0.4355  
 -0.3846 -0.0241  0.0785  0.0602 -1.0543  
  
i =  69LLF = -1283.6697step = 0.250Parameters  -3.4630  2.5308  1.0737 -1.9825 -2.1716 -1.0383 -2.1867  0.0201  0.1974  0.0763 -0.8583  0.0903  0.0055  0.4355  
 -0.3846 -0.0241  0.0785  0.0602 -1.0543  
  
i =  70LLF = -1283.6697step = 0.500Parameters  -3.4630  2.5308  1.0737 -1.9825 -2.1716 -1.0383 -2.1867  0.0201  0.1974  0.0763 -0.8583  0.0903  0.0055  0.4355  
 -0.3846 -0.0241  0.0785  0.0602 -1.0543  
  
i =  71LLF = -1283.6697step = 0.500Parameters  -3.4630  2.5308  1.0737 -1.9825 -2.1716 -1.0383 -2.1867  0.0201  0.1974  0.0763 -0.8583  0.0903  0.0055  0.4355  
 -0.3846 -0.0241  0.0785  0.0602 -1.0543  
  
i =  72LLF = -1283.6697step = 1.000Parameters  -3.4630  2.5308  1.0737 -1.9825 -2.1716 -1.0383 -2.1867  0.0201  0.1974  0.0763 -0.8583  0.0903  0.0055  0.4355  
 -0.3846 -0.0241  0.0785  0.0602 -1.0543  
  
Final results:
---------------------------------------------------------- 
 Variables      Value    Std.Err    Z-Value    P-Value  
---------------------------------------------------------- 
CONSTANT       -3.46304    0.98377   -3.52016    0.00043  
NEWENG          2.53079    0.68701    3.68375    0.00023  
MA              1.07372    0.34623    3.10114    0.00193  
ENC            -1.98249    0.63591   -3.11756    0.00182  
WNC            -2.17160    0.82481   -2.63284    0.00847  
PACIFIC        -1.03826    0.30497   -3.40444    0.00066  
Mount          -2.18668    0.65103   -3.35881    0.00078  
METRO           0.02014    0.16850    0.11952    0.90486  
MICRO           0.19736    0.20729    0.95208    0.34106  
HDD302          0.07633    0.08784    0.86895    0.38487  
CDD302         -0.85831    0.47164   -1.81983    0.06878  
HDDCDD          0.09032    0.11203    0.80624    0.42010  
HOUSEAGE        0.00548    0.00225    2.43530    0.01488  
OWN             0.43553    0.16937    2.57140    0.01013  
MOBILE         -0.38462    0.23248   -1.65441    0.09804  
INCOME_het     -0.02412    0.04252   -0.56732    0.57050  
HDD302_het      0.07850    0.03819    2.05564    0.03982  
HOUSEAGE_het    0.06019    0.01340    4.49060    0.00001  
NEWENG_het     -1.05435    0.45556   -2.31439    0.02065  
  

The likelihood ratio value is: 16.3314
This is the Prob. of Type I error:  0.05  
This is the critical Chi-Sq. Stat:   9.488-->Reject H0 

Question 2(b) - Compute average elasticity of HDD302 using the heteroskedastic model

diff = elast_hdd_het(b_phet);
fprintf('\nThe average elasticity of HDD302 is: %8.4f\n', diff);

% Delta Method to obtain variance of discrete change effect
diff_der = Grad(b_phet,'elast_hdd_het',1);
diff_cov = diff_der * c_phet * diff_der';
diff_se = sqrt(diff_cov);
fprintf('\nThe s.e. for the average elasticity of HDD302 is: %8.4f\n', diff_se);

% perform t-test (elasticity = 0)
tstat = (diff - 0)/sqrt(diff_cov);
pvalue = 2*(1-tcdf(abs(tstat),numobs-numc));
fprintf('\nThe t-stat for testing average elasticity of HDD302 = 0 is %8.4f (p-value = %5.4f)\n',...
        tstat, pvalue);
if pvalue<0.05
    disp('Reject H_0 under 0.05 significance level');
else
    disp('Fail to reject H_0 under 0.05 significance level');
end
The average elasticity of HDD302 is:   3.1722

The s.e. for the average elasticity of HDD302 is:   1.1433

The t-stat for testing average elasticity of HDD302 = 0 is   2.7746 (p-value = 0.0055)
Reject H_0 under 0.05 significance level

end of code

diary off;