setwd("d:/r") # ディレクトリーをd:/r にする

# 生産関数データの読み込み
product <- read.csv("product.csv", header=T) # データの頭にタイトルがあるので T とする
product # データを表示

# 各データの読み込みとデータ加工
Y <- product[,1] ; LY=log(Y) # データ Y は１列目にある
K <- product[,2] ; LK=log(K) # データ K は２列目にある
L <- product[,3] ; LL=log(L) # データ L は３列目にある

# ライブラリーの呼出
library(lmtest) # bptest、dwtest で用いる
library(tseries) # jarque.bera.test で用いる

# -------- log(Y) = a + b*log(K) + c*log(L) + u --------------
# 最小２乗法
result <- lm(LY~LK+LL) # LY=a+b*LK+c*LL+u
summary(result) # 結果を表示
deviance(result)  # 残差変動

# 誤差項の仮定をテストする統計量
fv <- result$fitted.values # 計算値ベクトルを fv とおく。
bptest(result,~fv**2) # LM 検定（分散均一性）
dwtest(result) # Durbin-Watson 検定（系列相関）
e <- residuals(result) # 残差ベクトルを e とおく。
jarque.bera.test(e) # Jarque-Bera 検定（正規性）
resettest(result,power=2:3) # Ramsey's RESET test （関数型の誤りがあるか）
# 上記の場合は、計算値の２乗と３乗を使って回帰している。もし２乗だけなら次式のようにする。
# resettest(result,power=2)
logLik(result) # 最大対数尤度
AIC(result) # 赤池情報量基準（定数項、説明変数２個、誤差項の分散を含め、４で計算）

plot(result)  # 判断のためのグラフがいくつか表示される。

# -------- log(Y/L) = a + b*log(K/L) + u --------------
# 一次同次関数の推定
LYL=log(Y/L)
LKL=log(K/L)
result <- lm(LYL~LKL)
summary(result) # 結果を表示
deviance(result)  # 残差変動

# 誤差項の仮定をテストする統計量
fv <- result$fitted.values # 計算値ベクトルを fv とおく。
bptest(result,~fv**2) # LM 検定（分散均一性）
dwtest(result) # Durbin-Watson 検定（系列相関）
e <- residuals(result) # 残差ベクトルを e とおく。
jarque.bera.test(e) # Jarque-Bera 検定（正規性）
reset(result) # Ramsey's RESET test （関数型の誤りがあるか）
logLik(result) # 最大対数尤度
AIC(result) # 赤池情報量基準（定数項、説明変数、誤差項の分散を含め、３で計算）