# Heteroskedasticity case
# Modeling food expenditure on disposible income

setwd("D:/Lectures/Econometrics/Heteroskedasticity/")

data	<- as.matrix(read.table(file="GHJ.txt"))

n	<- nrow(data)
p	<- ncol(data)-1

y	<- data[,1]					#vector of dependent var.
X	<- cbind(matrix(1,n,1),data[,2:(p+1)])	#matrix of explanatory var. and unit con.

# OLS for weight matrix
## residual
I	<- diag(1,n,n)				#Identity matrix
H	<- X%*%solve(t(X)%*%X)%*%t(X)		#Hat matrix

e	<- (I-H)%*%y				#Residual


## ploting residual vs disposible income
plot(data[,2], e, main="Plot residual melawan X", 
   xlab="Disposible income", ylab="Residual", pch=19)

## breuscch-pagan test
yh	<- e^2
e_h	<- (I-X%*%solve(t(X)%*%X)%*%t(X))%*%yh
R2_h	<- 1-(t(e_h)%*%e_h)/(n-1)/var(yh)
R2_h

ksi_h	<- n*R2_h
ksi_h
pv_h	<- 1-pchisq(ksi_h,p)
pv_h

# WLS
e2_l	<- as.vector(log(e^2))					#
yw<- e2_l
alpha<-as.matrix(solve(t(X)%*%X)%*%t(X)%*%yw)
alphaz<-X%*%alpha

w	<- as.vector(1/sqrt(exp(alphaz))) #1/h(i)
W	<- diag(w)					#psi^-1

b	<- solve(t(X)%*%X)%*%t(X)%*%y
b_w	<- solve(t(X)%*%W%*%X)%*%t(X)%*%W%*%y

e_w	<- sqrt(W)%*%(y-X%*%b_w)

s2	<- (t(e)%*%e)/(n-p-1)
s2_w	<- (t(e_w)%*%e_w)/(n-p-1)

se_b	<- sqrt(s2*diag(solve(t(X)%*%X)))
se_bw	<- sqrt(s2_w*diag(solve(t(X)%*%W%*%X)))

se_wh8i	<- sqrt(diag(solve(t(X)%*%X)%*%(t(X)%*%diag(as.vector(e^2))%*%X)%*%solve(t(X)%*%X)))