date: “Last compiled on 09 October, 2023”

2.1 几个例子

library(quantmod)
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: TTR
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
AAPL <- getSymbols("AAPL", src="tiingo", auto.assign=FALSE, api.key="9667b4325f683d9d0ae0734af175d886eceb7781")
chartSeries(
  AAPL$AAPL.Adjusted, type="line", TA=NULL,
  subset="2017/2023",
  major.ticks="years", minor.ticks=FALSE,
  theme="white", name="Apple"
  )

library(readr)
da <- read_table(
  "D:/齐安静 教学/时间序列分析/北大/ftsdata/q-ko-earns8309.txt",
  col_types=cols(
    .default = col_double(),
    pends = col_date("%Y%m%d"),
    anntime = col_date("%Y%m%d")
    ) )
ko.Rqtr <- xts(da[["value"]], da[["pends"]])
head(ko.Rqtr)
##              [,1]
## 1983-03-31 0.0375
## 1983-06-30 0.0492
## 1983-09-30 0.0463
## 1983-12-31 0.0379
## 1984-03-31 0.0425
## 1984-06-30 0.0583
chartSeries(
  ko.Rqtr, type="line", TA=NULL,
  major.ticks="years", minor.ticks=FALSE,
  theme="white", name="Coca Kola Quarterly Return"
  )

- 不同的季节用不同的颜色

library(lubridate)
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
tmp.x <- year(index(ko.Rqtr)) + (quarter(index(ko.Rqtr))-1)/4
tmp.y <- c(coredata(ko.Rqtr))
plot(tmp.x, tmp.y, type="l", col="gray",
     xlab="year", ylab="Return")
cpal <- c("green", "red", "yellow", "black")
points(tmp.x, tmp.y, pch=16, 
       col=cpal[quarter(index(ko.Rqtr))])
legend("topleft", pch=16, col=cpal,
       legend=c("Spring", "Summer", "Autumn", "Winter"))

- ggplot2绘制图形

library(magrittr)
library(dplyr)
## 
## ######################### Warning from 'xts' package ##########################
## #                                                                             #
## # The dplyr lag() function breaks how base R's lag() function is supposed to  #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
## # source() into this session won't work correctly.                            #
## #                                                                             #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
## # dplyr from breaking base R's lag() function.                                #
## #                                                                             #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## ###############################################################################
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:xts':
## 
##     first, last
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
da_b <- da %>%
  mutate(time = as.yearqtr(pends),
    qtr = factor(quarter(pends),
       levels = 1:4,
       labels=c("Spring", "Summer", "Autumn", "Winter")))
ggplot(data = da_b, mapping = aes(
  x = time, y = value)) +
  geom_point(mapping=aes(color=qtr)) +
  geom_line()

d <- read_table(
  "D:/齐安静 教学/时间序列分析/北大/ftsdata/m-ibmsp-2611.txt",
  col_types=cols(.default=col_double(),
                 date=col_date(format="%Y%m%d")))
## Warning: 1 parsing failure.
## row col  expected    actual                                                        file
## 226  -- 3 columns 4 columns 'D:/齐安静 教学/时间序列分析/北大/ftsdata/m-ibmsp-2611.txt'
head(d)
## # A tibble: 6 × 3
##   date           ibm       sp
##   <date>       <dbl>    <dbl>
## 1 1926-01-30 -0.0104  0.0225 
## 2 1926-02-27 -0.0245 -0.0440 
## 3 1926-03-31 -0.116  -0.0591 
## 4 1926-04-30  0.0898  0.0227 
## 5 1926-05-28  0.0369  0.00768
## 6 1926-06-30  0.0685  0.0432
sp.rmon <- xts(log(1 + d[["sp"]]), d[["date"]])
chartSeries(
  sp.rmon, type="line", TA=NULL,
  major.ticks="auto", minor.ticks=FALSE,
  theme="white", name="S&P 500 Monthly Returns"
  )

- 美国国债3月期和6月期周利率

library(readr)
d <- read_table(
  "D:/齐安静 教学/时间序列分析/北大/ftsdata/w-tb3ms.txt",
  col_types=cols(.default=col_double()))
library(lubridate)
library(xts)
x1 <- xts(d[["rate"]], make_date(d[["year"]], d[["mon"]], d[["day"]]))
d2 <- read_table(
  "D:/齐安静 教学/时间序列分析/北大/ftsdata/w-tb6ms.txt",
  col_types=cols(.default=col_double()))
x2 <- xts(d2[["rate"]], make_date(d2[["year"]], d2[["mon"]], d2[["day"]]))
tb36ms <- merge(x1, x2)
names(tb36ms) <- c("tb3ms", "tb6ms")
df <- data.frame(time=time(tb36ms), tb3ms=tb36ms[,1], tb6ms=tb36ms[,2])
# 绘制图形

plot(df$time, df$tb3ms, type='l',  
     col='black', xlab='', ylab='')

lines(df$time, df$tb6ms, col='red') 
grid()
legend("topright", legend = c("tb3ms","tb6ms"), col=c("black","red"), lty=1:2)

plot(tb36ms, type="l", subset="2004")

plot(tb36ms, type="l", subset="1980")

2.2 平稳性

da = read.table ("D:/齐安静 教学/时间序列分析/北大/ftsdata/m-ibmsp6709.txt", header=T)
head(da)
##       date       ibm        sp
## 1 19670131  0.075370  0.078178
## 2 19670228  0.079099  0.001963
## 3 19670331  0.048837  0.039410
## 4 19670428  0.100887  0.042239
## 5 19670531 -0.035234 -0.052441
## 6 19670630  0.067024  0.017512
ibm=da$ibm
sp5=da$sp
cor(sp5, ibm)
## [1] 0.5856544
cor(sp5 ,ibm, method= "spearman")
## [1] 0.5860817
cor(sp5, ibm, method="kendall")
## [1] 0.4196587
par(pty="s") # 把图片设定成正方形
plot(sp5,ibm, 
     xlab="S&P index returns", ylab="IBM stock returns",  xlim = c(-0.2,0.2))

### 2.2.1 自相关与白噪声 - 案例

d <- read_table(
  "D:/齐安静 教学/时间序列分析/北大/ftsdata/m-dec12910.txt", col_types=cols(
    .default=col_double(),
    date=col_date(format="%Y%m%d")))
head(d)
## # A tibble: 6 × 5
##   date           dec1    dec2     dec9    dec10
##   <date>        <dbl>   <dbl>    <dbl>    <dbl>
## 1 1967-01-31  0.0686   0.0804  0.181   0.212   
## 2 1967-02-28  0.00874  0.0110  0.0488  0.0649  
## 3 1967-03-31  0.0397   0.0354  0.0675  0.0689  
## 4 1967-04-28  0.0440   0.0375  0.0408  0.0446  
## 5 1967-05-31 -0.0506  -0.0362 -0.00219 0.000295
## 6 1967-06-30  0.0150   0.0189  0.102   0.119
library(xts)
dec <- xts(as.matrix(d[,-1]), d$date)
tclass(dec) <- "yearmon"
d10 <- ts(coredata(dec)[,"dec10"], start=c(1967,1), frequency=12)
plot(d10, main="CRSP Lower 10% Mothly Returns")

acf(d10)

d10=d$dec10 # select the Decile 10 returns
dec10=ts(d10,frequency=12,start=c(1967,1))
par(mfcol=c(1,1))
plot(dec10,xlab="year",ylab="returns")
title(main="(a): Simple returns")

acf(d10,lag=24) # command to obtain sample ACF of the data

TSA::acf(d10)

forecast::Acf(d10)
## Registered S3 methods overwritten by 'forecast':
##   method       from
##   fitted.Arima TSA 
##   plot.Arima   TSA

- 案例 （续）小市值股票倾向于在每年的一月份有正的收益率

tmp1 <- acf(d10, plot=FALSE)
tmp1$lag
## , , 1
## 
##       [,1]
##  [1,]    0
##  [2,]    1
##  [3,]    2
##  [4,]    3
##  [5,]    4
##  [6,]    5
##  [7,]    6
##  [8,]    7
##  [9,]    8
## [10,]    9
## [11,]   10
## [12,]   11
## [13,]   12
## [14,]   13
## [15,]   14
## [16,]   15
## [17,]   16
## [18,]   17
## [19,]   18
## [20,]   19
## [21,]   20
## [22,]   21
## [23,]   22
## [24,]   23
## [25,]   24
## [26,]   25
## [27,]   26
## [28,]   27
tmp1 <- acf(d10, plot=FALSE)
tmp1
## 
## Autocorrelations of series 'd10', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000  0.227 -0.019 -0.021  0.011  0.003 -0.028 -0.017 -0.049 -0.040  0.013 
##     11     12     13     14     15     16     17     18     19     20     21 
##  0.061  0.130 -0.037 -0.082 -0.021  0.017 -0.014 -0.059 -0.082 -0.064 -0.040 
##     22     23     24     25     26     27 
##  0.018 -0.015  0.052 -0.008  0.037  0.015
r12 <- tmp1$acf[abs(tmp1$lag-12)<1E-10] 
# 自动选择滞后12期
# 其中tmp1$lag，是滞后的期数，
# 减去12是一个新的数列
# 新的序列中，小于非常小的数1E-10，意味着就是要找到原来的滞后12阶才可以。
r12
## [1] 0.130411
t12 <- sqrt(tmp1$n.used)*r12; t12
## [1] 2.962369
pv <- 2*(1 - pnorm(abs(t12))); pv
## [1] 0.003052812
# 教材中的代码
f1=acf(d10,lag=24)

f1$acf
## , , 1
## 
##               [,1]
##  [1,]  1.000000000
##  [2,]  0.227386585
##  [3,] -0.019026447
##  [4,] -0.021258247
##  [5,]  0.011011345
##  [6,]  0.002676057
##  [7,] -0.027654887
##  [8,] -0.016910608
##  [9,] -0.049183690
## [10,] -0.039617756
## [11,]  0.013265549
## [12,]  0.061013220
## [13,]  0.130411045
## [14,] -0.036881195
## [15,] -0.082462743
## [16,] -0.020950139
## [17,]  0.016726386
## [18,] -0.013961209
## [19,] -0.059422809
## [20,] -0.082246074
## [21,] -0.063641596
## [22,] -0.039858376
## [23,]  0.017770989
## [24,] -0.015413528
## [25,]  0.052212082
tt=f1$acf[13]*sqrt(516)
tt
## [1] 2.962369

案例 检验IBM股票月收益率是否白噪声

library(readr)
d <- read_table(
  "D:/齐安静 教学/时间序列分析/北大/ftsdata/m-ibmsp-2611.txt", col_types=cols(
    .default=col_double(),
    date=col_date(format="%Y%m%d")))
## Warning: 1 parsing failure.
## row col  expected    actual                                                        file
## 226  -- 3 columns 4 columns 'D:/齐安静 教学/时间序列分析/北大/ftsdata/m-ibmsp-2611.txt'
ibm <- ts(d[["ibm"]], start=c(1926,1), frequency=12)
forecast::Acf(ibm)

- 简单收益率的白噪声检验

Box.test(ibm, lag=12, type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  ibm
## X-squared = 13.098, df = 12, p-value = 0.362
Box.test(ibm, lag=24, type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  ibm
## X-squared = 35.384, df = 24, p-value = 0.0629
  • 对数收益率的白噪声检验
Box.test(log(1 + ibm), lag=12, type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  log(1 + ibm)
## X-squared = 12.814, df = 12, p-value = 0.3827
Box.test(log(1 + ibm), lag=24, type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  log(1 + ibm)
## X-squared = 34.506, df = 24, p-value = 0.07607
  • 对10分位组合收益率的白噪声检验
forecast::Acf(d10)

Box.test(d10, type="Ljung-Box", lag=12)
## 
##  Box-Ljung test
## 
## data:  d10
## X-squared = 41.06, df = 12, p-value = 4.789e-05
Box.test(d10, type="Ljung-Box", lag=24)
## 
##  Box-Ljung test
## 
## data:  d10
## X-squared = 56.246, df = 24, p-value = 0.0002122