?>

 

#  read csv

A  <- read.csv(file=’Stock.csv’)

str(A)

head(A)

 

windows()

plot(A)

 

# convert to time series

B <- ts(A,  frequency = 1)     ###ts時間數列

str(B)     ###B是Time Series

windows()

plot(B)

 

# extract SIC2454

price <- B[,5]

windows()

plot(price)

###在左框輸入str(price)看其結構

 

return <- diff(price)

windows()

plot(return)

 

rate <- diff(log(price))

windows()

plot(rate)

 

#install.packages(‘forecast’)

# automatically check ARMIA

library(forecast)

 

auto.arima(price)

auto.arima(return)

auto.arima(rate)

 

#install.packages(‘tseries’)

# unit root check

library(tseries)

adf.test(price,k=0)

adf.test(return,k=0)

adf.test(rate,k=0)

 

###############  ACF PACF

windows()

acf(price)

windows()

pacf(price)

 

 

windows()

acf(return)

windows()

pacf(return)

 

 

windows()

acf(rate)

windows()

pacf(rate)

 

 

 

#### plot stock price

windows()

autoplot(price)

windows()

autoplot(return)

windows()

autoplot(rate)

 

 

########## redirect + acf+pacf

windows()

price %>% diff() %>% ggtsdisplay(main=”Return”)

 

windows()

rate %>%  ggtsdisplay(main=”Rate of Return”)      ###價格做差分就是報酬  ###價格做差分以後的數列就是定態

 

#### model

fit1 <- Arima(price, order=c(0,1,0))

fit1

checkresiduals(fit1)

auto.arima(fit1$residuals)

 

 

windows()

hist(fit1$residuals,breaks=30)

 

windows()

library(“car”)

qqPlot(fit1$residuals)

Box.test(fit1$residuals,type=’Ljung’,lag=1)

Box.test(fit1$residuals,type=’Ljung’,lag=6)

 

##### forecast

windows()

autoplot(forecast(fit1))

forecast(fit1)

 

 

 

######## another case

 

auto.arima(rate)

 

fit2 <- Arima(rate, order=c(1,1,1))

fit2

 

windows()

checkresiduals(fit2)

auto.arima(fit2$residuals)

windows()

acf(fit2$residuals)

windows()

pacf(fit2$residuals)

 

windows()

qqPlot(fit2$residuals)

Box.test(fit2$residuals,type=’Ljung’,lag=6)

windows()

autoplot(forecast(fit2),main=’XXX’)

 

 

 

 

221208 計量經濟學
?>