Econometric analysis Whilliam H.Greene Prentice hall

```{r}

#γ分布の最尤推定 p.71

#勾配上昇法

XX=c(3,1,4,5,3)

n=length(XX)

X=rep(1,2)

eta=10^(-2)

ite=4*10^4

for(l in 1:ite){

X1=-n*digamma(X[1])-n*log(X[2])+(X[1]-1)*sum(log(XX))

X2=-n*X[1]/X[2]+sum(XX)/(X[2]^2)

X[1]=X[1]+eta*X1

X[2]=X[2]+eta*X2

loglik=-n*log(gamma(X[1]))-n*X[1]*log(X[2])+(X[1]-1)*mean(log(XX))-sum(XX)/X[2]

#print(loglik)

}


```



```{r}

#γ分布の最尤推定 p.71

XX=c(3,1,4,5,3)

n=length(XX)


#fiacco,mccormick(1968),murtagh,sargent(1969)

f=function(x){

f1=-n*digamma(x[1])-n*log(x[2])+(x[1]-1)*sum(log(XX))

f2=-n*x[1]/x[2]+sum(XX)/(x[2]^2)

return(c(f1,f2))

}

X=rep(1,2)

H=diag(1,length(X))

ite=3000

eta=10^(-2)

h=0.01

for(l in 1:ite){

f_val=f(X)

sigma=t(t(-H%*%f_val))

X=X+eta*sigma

y=f(X)-f_val

H=H+(eta*sigma-H%*%y)%*%t(eta*sigma-H%*%y)/as.numeric(t(y)%*%(eta*sigma-H%*%y))

loglik=-n*log(gamma(X[1]))-n*X[1]*log(X[2])+(X[1]-1)*mean(log(XX))-sum(XX)/X[2]

#print(f(X))

print(loglik)

}







```


```{r}

#γ分布の最尤推定　p.71

XX=c(3,1,4,5,3)

n=length(XX)


#Broyden Quasi-Newton Algorithm 

f=function(x){

f1=-n*digamma(x[1])-n*log(x[2])+(x[1]-1)*sum(log(XX))

f2=-n*x[1]/x[2]+sum(XX)/(x[2]^2)

return(c(f1,f2))

}

X=rep(1,2)

ite=2000

eta=0.01

H=diag(rep(0.1,length(X)))

for(l in 1:ite){

X_pre=X  

X=X-eta*H%*%f(X)  

s=X-X_pre

y=f(X)-f(X_pre)

H=H+((s-H%*%y)/as.numeric(t(s)%*%H%*%y))%*%t(s)%*%H

#print((f(X)))

loglik=-n*log(gamma(X[1]))-n*X[1]*log(X[2])+(X[1]-1)*mean(log(XX))-sum(XX)/X[2]

print(loglik)


}  




```

```{r}

#Econometric analysis 

#p.304

data=data.frame(num=1:20,y=c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56))

Y=data$y

X=cbind(rep(1,nrow(X)),as.matrix(data[,colnames(data) %in% c("x1","x2","x3")]))

n=nrow(X)

sig=1

beta=rep(1,ncol(X))

lam=1

ite=10^5

eta=10^(-5)

for(l in 1:ite){
  
beta=beta+eta*(t(X)%*%(Y-X%*%beta)/as.numeric(sig)-2*beta)  
  
sig=sig+eta*(-n/(2*sig)+t(Y-X%*%beta)%*%(Y-X%*%beta)/as.numeric(2*sig^2))

lam=lam+eta*(-sum(beta^2))

loglik=-n*(log(2*pi)+log(sig))/2-t(Y-X%*%beta)%*%(Y-X%*%beta)/as.numeric(2*sig^2)-lam*2*sum(beta^2)
 
print(as.numeric(loglik))
 
}

```

```{r}


#Econometric analysis 

#p.474

data=data.frame(num=1:20,y=c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56))

Y=data$y

X=cbind(rep(1,nrow(X)),as.matrix(data[,colnames(data) %in% c("x1","x2","x3")]))

data(anscombe)

Y=anscombe$y4*100

X=cbind(rep(1,length(Y)),anscombe[,colnames(anscombe) %in% c("y1","y2","y3")])

X=as.matrix(X)

n=nrow(X)

sig=1

beta=rep(1,ncol(X))

sita=1

ite=10^6

eta=10^(-2)

for(l in 1:ite){

ep=log(Y)+sita*Y-X%*%beta

dbeta=t(X)%*%ep/sig

beta=beta+eta*dbeta

dsita=sum(Y/(1+sita*Y)-Y*ep/sig)

dsig=sum((ep^2/sig-1)/(2*sig))

sita=sita+eta*dsita

sig=sig+eta*dsig

loglik=sum(log(1+sita*Y)-log(Y)-(log(2*pi)+log(sig))/2-(log(Y)+sita*Y-X%*%beta)/(2*sig))

print(as.numeric(loglik))

}

plot(c(1:length(Y)),log(Y)+sita*Y,xlim=c(1,length(Y)),ylim=c(min(log(Y)+sita*Y),max(log(Y)+sita*Y)),col=2,xlab="samples",ylab="values")

par(new=T)

plot(c(1:length(Y)),(X)%*%beta,xlim=c(1,length(Y)),ylim=c(min(log(Y)+sita*Y),max(log(Y)+sita*Y)),col=3,xlab="samples",ylab="values")




```


```{r}

#p.476のデータ

Y=c(126.148,3201.486,690.67,56.296,304.531,723.028,992.169,35.796,494.515,124.948,73.328,29.467,415.262,241.53,4079.554,652.085,667.113,940.43,1611.899,617.579,527.413,174.394,636.948,22.7,349.711)

X=matrix(c(c(3.804,185.446,39.712,6.547,11.53,58.987,112.884,2.698,10.36,5.213,3.763,1.967,17.546,15.347,435.105,32.84,33.292,72.974,157.978,34.324,22.736,7.173,30.807,1.543,22.001),c(31.551,452.844,124.074,19.181,45.534,88.391,148.53,8.017,86.189,12,15.9,6.47,69.342,39.416,490.384,84.831,83.033,190.094,259.916,98.152,109.728,31.301,87.963,4.063,52.818)),ncol=2)


#Econometric analysis 

#p.484のモデル BOX-COX trans


n=nrow(X)

sig=1

beta=rep(0.001,ncol(X))

lam=0.9

ite=10^6

eta=10^(-8)

for(l in 1:ite){

Y_lam=(Y^lam-1)/lam

X_lam=(X^lam-1)/lam
  
ep=Y_lam-(X_lam)%*%beta

dbeta=t(X_lam)%*%ep/sig

beta=beta+eta*dbeta

dlam=sum(log(Y)-ep*(((Y^lam)*log(Y)-Y_lam)/lam-(1/lam)*((X^lam)*log(X)-X_lam)%*%beta)/sig)

dsig=sum((ep^2/sig-1)/(2*sig))

lam=lam+eta*dlam

sig=sig+eta*dsig

loglik=(lam-1)*sum(log(Y))-log(2*pi)/2-length(Y)*log(sig)/2-sum(ep^2)/(2*sig)

#print(as.numeric(loglik))

}

plot(c(1:length(Y)),Y_lam,xlim=c(1,length(Y)),ylim=c(min(Y_lam),max(Y_lam)),col=2,xlab="samples",ylab="values")

par(new=T)

plot(c(1:length(Y)),(X_lam)%*%beta,xlim=c(1,length(Y)),ylim=c(min(Y_lam),max(Y_lam)),col=3,xlab="samples",ylab="values")



```

```{r}

#p.476のデータ

Y=c(126.148,3201.486,690.67,56.296,304.531,723.028,992.169,35.796,494.515,124.948,73.328,29.467,415.262,241.53,4079.554,652.085,667.113,940.43,1611.899,617.579,527.413,174.394,636.948,22.7,349.711)

X=matrix(c(rep(1,25),c(3.804,185.446,39.712,6.547,11.53,58.987,112.884,2.698,10.36,5.213,3.763,1.967,17.546,15.347,435.105,32.84,33.292,72.974,157.978,34.324,22.736,7.173,30.807,1.543,22.001),c(31.551,452.844,124.074,19.181,45.534,88.391,148.53,8.017,86.189,12,15.9,6.47,69.342,39.416,490.384,84.831,83.033,190.094,259.916,98.152,109.728,31.301,87.963,4.063,52.818)),ncol=3)


#Econometric analysis 

#p.484のモデル


n=nrow(X)

sig=1

beta=rep(0.001,ncol(X))

lam=0.9

ite=10^6

eta=10^(-7)

for(l in 1:ite){

ep=Y^lam-(X^lam)%*%beta

dbeta=t(X^lam)%*%ep/sig

beta=beta+eta*dbeta

dlam=sum(log(Y)-ep*((Y^lam)*log(Y)-((X^lam)*log(X))%*%beta)/sig)

dsig=sum((ep^2/sig-1)/(2*sig))

lam=lam+eta*dlam

sig=sig+eta*dsig

loglik=(lam-1)*sum(log(Y))-log(2*pi)/2-length(Y)*log(sig)/2-sum(ep^2)/(2*sig)

#print(as.numeric(loglik))

}

plot(c(1:length(Y)),Y^lam,xlim=c(1,length(Y)),ylim=c(min(Y^lam),max(Y^lam)),col=2,xlab="samples",ylab="values")

par(new=T)

plot(c(1:length(Y)),(X^lam)%*%beta,xlim=c(1,length(Y)),ylim=c(min(Y^lam),max(Y^lam)),col=3,xlab="samples",ylab="values")



```

```{r}



#Econometric analysis 

#p.484

data=data.frame(num=1:20,y=c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56))

Y=data$y

X=cbind(rep(1,nrow(X)),as.matrix(data[,colnames(data) %in% c("x1","x2","x3")]))

data(anscombe)

Y=anscombe$y4*100

X=cbind(rep(1,length(Y)),anscombe[,colnames(anscombe) %in% c("y1","y2","y3")])

X=as.matrix(X)

n=nrow(X)

sig=1

beta=rep(1,ncol(X))

lam=0.9

ite=10^6

eta=10^(-4)

for(l in 1:ite){

ep=Y^lam-(X^lam)%*%beta

dbeta=t(X^lam)%*%ep/sig

beta=beta+eta*dbeta

dlam=sum(log(Y)-ep*((Y^lam)*log(Y)-((X^lam)*log(X))%*%beta)/sig)

dsig=sum((ep^2/sig-1)/(2*sig))

lam=lam+eta*dlam

sig=sig+eta*dsig

loglik=(lam-1)*sum(log(Y))-log(2*pi)/2-length(Y)*log(sig)/2-sum(ep^2)/(2*sig)

print(as.numeric(loglik))

}

plot(c(1:length(Y)),Y^lam,xlim=c(1,length(Y)),ylim=c(min(Y^lam),max(Y^lam)),col=2,xlab="samples",ylab="values")

par(new=T)

plot(c(1:length(Y)),(X^lam)%*%beta,xlim=c(1,length(Y)),ylim=c(min(Y^lam),max(Y^lam)),col=3,xlab="samples",ylab="values")



```

```{r}

#Econometric analysis 

#p.513 FGLS estimator

data=data.frame(num=1:20,y=c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56))

Y=data$y

X=cbind(rep(1,nrow(X)),as.matrix(data[,colnames(data) %in% c("x1","x2","x3")]))

n=nrow(X)

sig=1

beta=rep(1,ncol(X))

gam=diag(rep(1,nrow(X)))

ite=10^5

eta=10^(-7)

for(l in 1:ite){
  
beta=beta+eta*t(X)%*%gam%*%(Y-X%*%beta)/as.numeric(sig)  
  
sig=sig+eta*(-n/(2*sig)+t(Y-X%*%beta)%*%gam%*%(Y-X%*%beta)/as.numeric(2*sig^2))

gam=gam+eta*(solve(gam)-(Y-X%*%beta)%*%t(Y-X%*%beta)/as.numeric(sig))

loglik=-n*(log(2*pi)+log(sig))/2-t(Y-X%*%beta)%*%gam%*%(Y-X%*%beta)/as.numeric(2*sig^2)+log(det(gam))/2
 
print(as.numeric(loglik))
 
}

```


```{r}

#p.565 Multiplicative Heteroscedasticity

y=c(126.148,3201.486,690.67,56.296,304.531,723.028,992.169,35.796,494.515,124.948,73.328,29.467,415.262,241.53,4079.554,652.085,667.113,940.43,1611.899,617.579,527.413,174.394,636.948,22.7,349.711)

X=matrix(c(rep(1,25),c(3.804,185.446,39.712,6.547,11.53,58.987,112.884,2.698,10.36,5.213,3.763,1.967,17.546,15.347,435.105,32.84,33.292,72.974,157.978,34.324,22.736,7.173,30.807,1.543,22.001),c(31.551,452.844,124.074,19.181,45.534,88.391,148.53,8.017,86.189,12,15.9,6.47,69.342,39.416,490.384,84.831,83.033,190.094,259.916,98.152,109.728,31.301,87.963,4.063,52.818)),ncol=3)

Z=X[,2:3]

n=nrow(X)

ite=10^5

eta=10^(-3)

beta=rep(0.1,ncol(X))

alpha=rep(0.1,ncol(Z))

n=nrow(X)

loglik_vec=c()

for(l in 1:ite){
  
ep=y-X%*%beta  

inv_omega=diag(c(1/exp(Z%*%alpha)))
  
dbeta=t(X)%*%inv_omega%*%ep

dalpha=t(Z)%*%(ep^2/exp(Z%*%alpha)-1)/2

ddbeta=t(X)%*%inv_omega%*%X

ddalpha=(t(Z)%*%Z)/2

beta=beta+eta*solve(ddbeta)%*%dbeta

alpha=alpha+eta*solve(ddalpha)%*%alpha

ep=y-X%*%beta  

loglik=-n*log(2*pi)/2-sum(Z%*%alpha)/2-sum(ep^2/exp(Z%*%alpha))/2

#print(loglik)

loglik_vec=c(loglik_vec,loglik)
  
}


```

```{r}

#p.889

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22),sign=0)

Y=ifelse(data$y>0.5,1,0)

X=as.matrix(cbind(rep(1,nrow(data)),data[,colnames(data) %in% c("x1","x2","x3")]))

z1=c(1,0,1,0,1,1,0,0,1,1,0,0,1,0,1,1,0,1,1,1)

z2=abs(z1-1);z=cbind(z1,z2)

r=rep(0.01,ncol(z))

beta=rep(0.0001,ncol(X))

ite=10^5

eta=0.00001

for(l in 1:ite){
  
f=dnorm(X%*%beta/exp(z%*%r))  
  
fd=pnorm(X%*%beta/exp(z%*%r))

fd=ifelse(abs(fd)<10^(-15),10^(-15),fd)

fd=ifelse(abs(1-fd)<10^(-15),1-10^(-15),fd)

dbeta=t(X)%*%(exp(-z%*%r)*f*(Y-fd)/(fd*(1-fd)))
  
dr=t(c((-(X%*%beta)*exp(-z%*%r)*f*(Y-fd)/(fd*(1-fd)))))%*%z

beta=c(beta+eta*dbeta)

r=c(r+eta*dr)

loglik=sum(Y*log(fd)+(1-Y)*log(1-fd))

print(loglik)

}


```

```{r}

#p.919

X=as.matrix(iris[,!(colnames(iris) %in% "Species")])

signal=iris$Species

Species=unique(signal)

n=50

J=3

vec1=sample(1:10,n,replace=T)
vec2=sample(20:30,n,replace=T)
vec3=sample(40:60,n,replace=T)

d=c(vec1,vec2,vec3)

beta=rep(0.01,ncol(X))

ite=10^(4)

eta=10^(-5)

for(l in 1:ite){
  
mat=array(0,dim=c(n,J))  
  
for(s in 1:J){

mat[,s]=X[signal==Species[s],]%*%beta  
  
}

prob=exp(mat)/apply(exp(mat),1,sum)

ave=array(0,dim=c(n,ncol(X)))

for(s in 1:J){
  
ave=ave+X[signal==Species[s],]*c(prob[,s])
  
}

dbeta=array(0,dim=c(n,ncol(X)))

for(s in 1:J){
  
dbeta=dbeta+(X[signal==Species[s],]-ave)*c(d[signal==Species[s]])
  
}

dbeta=apply(dbeta,2,sum)

beta=beta+eta*dbeta

loglik=sum(matrix(d,ncol=J)*log(prob))

print(loglik)

}


d_mat=matrix(d,ncol=J)

propotion=t(t(d_mat)/apply(d_mat,1,sum))

plot(c(1:(n*J)),c(prob),xlim=c(1,n*J),ylim=c(0,1),col=2,xlab="index",ylab="prob")

par(new=T)

plot(c(1:(n*J)),c(propotion),xlim=c(1,n*J),ylim=c(0,1),col=3,xlab="index",ylab="prob")


```