#### 一：時間序列圖
# 商品房銷售價格指數的時間序列圖
library(psych)
#繪制時序圖
x<-ts(y, frequency =1, start=c(2004),end=c(2016))
plot(x)
#1階差分，並繪制出差分后序列的時序圖
y.dif<-diff(x)
plot(y.dif)
# 差分序列adf單位根檢驗
library(urca)
a<-ur.df(y.dif)
b<-summary(a)
taus<-b@cval[1,]
ts<-b@teststat[1:1]
# acf pacf
#繪制序列自相關圖和偏自相關圖
a<-acf(y.dif,plot=T)
a
b<-pacf(y.dif,plot=T)
b
# 采用AR(1)模型
y_l=y.dif[2:12]
x_l=y.dif[1:11]
ols_fit<-lm(y_l~x_l)
summary(ols_fit)
#ols下的殘差值
residuals(ols_fit)
#ols下的估計值
fitted(ols_fit)
#擬合曲線
plot(y_l~x_l)
abline(ols_fit)
### 局部回歸模型：
model=loess(y_l~x_l,control = loess.control(surface = "direct"),degree=2)
predictions1<- predict (model,x_l)
# 繪圖
plot(x_l,y_l,pch=19)
x1<-order(x_l)
nr<-length(x1)
x_arr<-array(1:nr)
y_arr<-array(1:nr)
for(iin c(1:nr)){
temp=x1[i]
x_arr[i]=x_l[temp]
y_arr[i]=predictions1[temp]
}
abline(ols_fit,col='red')
lines(x_arr,y_arr,col='blue')
#### 二： OLS 回歸
ols_fit<-lm(y~x)
summary(ols_fit)
#ols下的殘差值
residuals(ols_fit)
#ols下的預測值
fitted(ols_fit)
#擬合曲線
plot(y~x)
abline(ols_fit)
#### 三 NW核回歸
# 采用的是高斯核（Gaussian Kernal），故核函數是正態分布下的密度函數
kernalGaussian <- function(xData)
{
stdX <- sd(xData)
# 高斯寬帶的選擇：每個點處的最優帶寬
h <- 1.06*stdX*length(xData)^(-1/5)
print(h)
# 每個點處的核函數
kernalX <- 1/(h*sqrt(2*pi)) * exp(-xData^2/(2*h^2))
return(kernalX)
}
# Nadaraya-Waston核估計，參數xData , yData必須是矩陣，且長度一樣
kernalRegress <- function(xData , yData)
{
# 最終返回針對y的核回歸擬合的值
nData<-nrow(xData)
yRegress <- matrix(NaN , nrow = nData , ncol = 1)
for (iin c(1:nData))
{
x <- xData[i]
xXt <- matrix(x , nrow = nData, ncol = 1) - xData
# khx也就是權重
khX <- kernalGaussian(xXt)
# yRegress 加權算術平均值:求出x處的平均值
yRegress[i] <- sum(yData*khX)/sum(khX)
}
return(yRegress)
}
# 核回歸的檢測
#x,y排序
x1<-order(x)
nr<-length(x)
x_arr<-array(1:10)
y_arr<-array(1:10)
for(iin c(1:nr)){
temp=x1[i]
x_arr[i]=x[temp]
y_arr[i]=y[temp]
}
# 把x,y變成矩陣
x_matrix<- as.matrix(x_arr)
y_matrix<- as.matrix(y_arr)
# 核回歸
y_regress<-kernalRegress(x_matrix,y_matrix)
# 真實值和預測值
cbind(y_matrix,y_regress)
# 畫圖
plot(x_arr,y_matrix,xlab = "全體居民消費指數", ylab = "商品房銷售價格指數")
lines(x_arr,y_regress,col = 'red')
# 合並圖
plot(x_arr, y_arr, xlab = "全體居民消費指數", ylab = "商品房銷售價格指數", col = 1)
lines(x_arr, abline(ols_fit), lty = 1, col = 1)
lines(x_arr,y_regress, lty = 2, col = 2)
letters <- c("OLS model", "NW method")
legend("bottomright", legend = letters, lty = 1:2, col = 1:2, cex = 0.5)