setwd("d:/r")

# 生鮮魚介データ（時系列）の読み込み
fish <- read.csv("fish.csv", header=T)
fish

# 各データの読み込み
F1 <- fish[,2] # 魚介金額
F2 <- fish[,3] # 魚介数量
F3 <- fish[,4] # 魚介価格
F4 <- fish[,5] # 魚介物価指数
H  <- fish[,6] # 世帯人員 
C  <- fish[,7] # 消費支出
P  <- fish[,8] # 総合物価指数
 
# データ加工（対数）
LQ=log(F2)        # 数量
LY=log(C/P*100)   # 実質消費支出
LP=log(F3/F4*100) # 実質価格
LH=log(H)         # 世帯人員

# 最小２乗法
result <- lm(LQ~LY+LP+LH)
summary(result) # 結果を表示

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

# 誤差項の仮定をテストする統計量
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 （関数型の誤りがあるか）
AIC(result) # 赤池情報量基準

# 誤差項にＡＲ１を仮定してρの推定値を計算
rhodata<-arima (residuals (result), order = c(1,0,0), include.mean = FALSE) # 誤差項にAR(1)
rho<-coef(rhodata)
rho

# ρの推定値を用いて変数を変換
LQ1=LQ-rho*LQ
LY1=LY-rho*LY
LP1=LP-rho*LP
LH1=LH-rho*LH

# 変換された変数で最小２乗法
result1 = lm(LQ1~LY1+LP1+LH1)
summary(result1) # 結果を表示

# 定数項を変換し直す
crho<-coef(result1)
c<-crho[1]/(1-rho)
c