chapter1
qaj
2023-08-14
date: “Last compiled on 27 September, 2023”
1.1 资产收益率
1.1.1 简单收益率
- 苹果公司2011-12-02到2011-12-09的每日股票收盘价
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.2 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsd.apple <- tibble(
date=c(ymd("2011-12-02"),
ymd("2011-12-05") + days(0:4)),
price=c(389.70, 393.01, 390.95, 389.09, 390.66, 393.62)
)
knitr::kable(d.apple)| date | price |
|---|---|
| 2011-12-02 | 389.70 |
| 2011-12-05 | 393.01 |
| 2011-12-06 | 390.95 |
| 2011-12-07 | 389.09 |
| 2011-12-08 | 390.66 |
| 2011-12-09 | 393.62 |
- 计算单期简单收益率的函数：
simple.return <- function(x){
x <- as.vector(x)
c(NA, diff(x) / x[1:(length(x)-1)])
}- 计算简单收益率(时间单位是交易日，无交易的日期忽略)， 精确到百分数的小数点后2位。
d.apple <- d.apple %>%
mutate(SR1=round(100*simple.return(price), 2))
knitr::kable(d.apple)| date | price | SR1 |
|---|---|---|
| 2011-12-02 | 389.70 | NA |
| 2011-12-05 | 393.01 | 0.85 |
| 2011-12-06 | 390.95 | -0.52 |
| 2011-12-07 | 389.09 | -0.48 |
| 2011-12-08 | 390.66 | 0.40 |
| 2011-12-09 | 393.62 | 0.76 |
- 从2011-12-02到2011-12-09是周五到周五， 这一周时间的多期简单收益率为：
(d.apple[6,"price"] - d.apple[1,"price"]) /
d.apple[1,"price"]## price
## 1 0.010059021.1.2 连续复利收益率
- \[R' = (1 + \frac{R}{n})^n - 1\]
- 作业1：借吧的例子
- 1000元，每天0.2元利息，实际年利率是多少?官方声称是7.3%
r <- 100*((1+0.2/1000)^365-1)
r## [1] 7.572269- 取年连续复利\(r=0.10\),\(m=1,2,3,12,52,365,\infty\) ， 初始资产1元，多次付息并计入账户， 一年后净值分别为
library(tidyverse)
m <- c(1, 2, 3, 12, 52, 365)
mlab <- c(paste(m), "Inf")
y <- c((1 + 0.20/m)^m, exp(0.20))
knitr::kable(tibble(`付息次数`=mlab, `年末净值`=y))| 付息次数 | 年末净值 |
|---|---|
| 1 | 1.200000 |
| 2 | 1.210000 |
| 3 | 1.213630 |
| 12 | 1.219391 |
| 52 | 1.220934 |
| 365 | 1.221336 |
| Inf | 1.221403 |
- 给定一个价格序列，计算对数收益率序列的函数如下：
log.return <- function(x){
c(NA, diff(log(x)))
}- 苹果公司股票的日对数收益率：
d.apple <- d.apple %>%
mutate(LR1=round(100*log.return(price), 2))
knitr::kable(d.apple)| date | price | SR1 | LR1 |
|---|---|---|---|
| 2011-12-02 | 389.70 | NA | NA |
| 2011-12-05 | 393.01 | 0.85 | 0.85 |
| 2011-12-06 | 390.95 | -0.52 | -0.53 |
| 2011-12-07 | 389.09 | -0.48 | -0.48 |
| 2011-12-08 | 390.66 | 0.40 | 0.40 |
| 2011-12-09 | 393.62 | 0.76 | 0.75 |
1.1.3 现值分析
1.1.3.1 贷款还贷分析
例如，设\(L=100\)万元，\(n=120\) 个月(10年)，\(R=0.005\) ， 则：
loan_pm <- function(L = 100, n = 120, R = 0.005){
L * R / (1 - (1 + R)^(-n))
}
loan_pm()## [1] 1.110205loan_pm() * 120 - 100## [1] 33.2246((loan_pm() * 120) / 100)^(1/10) - 1## [1] 0.02910205(1+0.005)^12-1## [1] 0.06167781- 例如，设\(L=100\)万元， \(n=120\)个月(10年)，\(R=0.005\) ， 则：
loan_demo <- function(L = 100, n = 120, R = 0.005){
A <- L * R / (1 - (1 + R)^(-n))
alpha <- 1 + R
jvec <- 1:n
Bvec <- L * alpha^jvec * (alpha^(n - jvec) - 1) / (alpha^n - 1)
Ivec <- L * (alpha - 1) * (alpha^(n) - alpha^(jvec-1)) / (alpha^n - 1)
Pvec <- A - Ivec
data.frame(
"月份" = jvec,
"还款" = A,
"还本" = Pvec,
"利息" = Ivec,
"剩余" = Bvec
)
}
tmp.d <- loan_demo()
knitr::kable(tmp.d[seq(12, 120, by=12),], digits=2, row.names=FALSE)| 月份 | 还款 | 还本 | 利息 | 剩余 |
|---|---|---|---|---|
| 12 | 1.11 | 0.64 | 0.47 | 92.47 |
| 24 | 1.11 | 0.68 | 0.43 | 84.48 |
| 36 | 1.11 | 0.73 | 0.38 | 76.00 |
| 48 | 1.11 | 0.77 | 0.34 | 66.99 |
| 60 | 1.11 | 0.82 | 0.29 | 57.43 |
| 72 | 1.11 | 0.87 | 0.24 | 47.27 |
| 84 | 1.11 | 0.92 | 0.19 | 36.49 |
| 96 | 1.11 | 0.98 | 0.13 | 25.05 |
| 108 | 1.11 | 1.04 | 0.07 | 12.90 |
| 120 | 1.11 | 1.10 | 0.01 | 0.00 |
可以看出本金的偿还额越到后期越多， 而还款前期主要支付的是利息。
作业2：汽车贷款利率
f <- function(x) 3111.11*((1+x)^36-1) - (1+x)^36*100000*(x)
r <- uniroot(f,c(0.001,0.9999))$root
cat((1+r)^12-1)## 0.07835412- 作业3：上海买房问题
- 2019年6月，A先生在上海买了一套住房，贷款300万，贷款年利率是4.6%。贷款期限是20年，等额本息每月还款，当年7月第1次还款。现在是2023年9月，由于利率下行，A先生的贷款利率下降到了4.2%。请问，A先生从下一个月开始，每月的还款额减少了多少？_______(答案带单位“元”，例如：325.12元。)
- 第1步 求出月利率
fleft1 <- function(R) (1+R)^12 - 1.046
fleft2 <- function(x) (1+x)^12 - 1.042
mrate1 <- uniroot(fleft1, c(0.001, 0.99))$root
mrate2 <- uniroot(fleft2, c(0.001, 0.99))$root- 第2步 利用上面的函数来计算每月还款额。
tmp1 <- loan_demo(3000000,240,mrate1)
principal <- tmp1[51,5]
anuity <- tmp1[51,2]- 第3步，按照本金和新的利率计算每月的还款额
tmp2 <- loan_demo(principal,240-51,mrate2)
new_anuity <- tmp2[1,2]
answer <- anuity-new_anuity
cat("每月的还款额减少了",sprintf("%.2f", answer),"元")## 每月的还款额减少了 505.10 元1.1.3.2 回报率
例题
- 投资100元， 分5年的回报分别为\(20,22,25,30,40\)元。 求回报率。
解答
- R函数uniroot()可以用来求解一元连续函数的单个根。
x <- c(20, 22, 25, 30, 40)
fleft <- function(R) -100 + sum(x / (1+R)^(1:5))
uniroot(fleft, c(-0.99, 10.0), extendInt="downX")## $root
## [1] 0.1016542
##
## $f.root
## [1] -0.001301806
##
## $iter
## [1] 12
##
## $init.it
## [1] NA
##
## $estim.prec
## [1] 6.103516e-05例题
- 设某民间借贷10000元，借期1年， 采用“三分利”， 指月利率\(3\%\)， 需要每月付息，不计复利， 到期偿还本金与最后一个月的利息。 大约相当于年利率\(36\%\)， 用现值分析方法计算实际年利率。
解答
x <- c(rep(300, 11), 10300)
fleft <- function(R) -10000 + sum(x / (1+R/12)^(1:12))
uniroot(fleft, c(-0.99, 10.0), extendInt="downX")## $root
## [1] 0.3600024
##
## $f.root
## [1] -0.02024998
##
## $iter
## [1] 9
##
## $init.it
## [1] NA
##
## $estim.prec
## [1] 6.103516e-051.1.7 小节
exp(0.0446) - 1## [1] 0.04560953- 简单收益率
d.ibm <- read_table(
"D:/齐安静 教学/时间序列分析/北大/ftsdata/d-ibm-0110.txt",
col_types=cols(date=col_date(format="%Y%m%d"),
return=col_double()))
with(d.ibm, plot(
date, return, type="l",
xlab="Date", ylab=expression(R[t])))
- 对数收益率
with(d.ibm, plot(
date, log(1 + return), type="l",
xlab="Date", ylab=expression(r[t])))
1.2.2 到期收益率
f <- function(y) 2.5/y*(1 - (1+y)^-6) + 100*(1+y)^-6
curve(f(x), 0.028, 0.052, xlab="半年期到期收益率", ylab="售价")
abline(v=c(0.030, 0.035, 0.040, 0.045, 0.050), col="gray")
pr <- c(97.29, 94.67, 92.14, 89.68, 87.31)
abline(h=pr, col="gray")
- P6中的表格。
library(tidyverse)
YTM <- c(6:10)
m <- 2
y <- YTM/(100*m)
F=100; alpha=0.05; n=3
k<-m*n
p <- F*alpha/(m*y)*(1 - (1+y)^(-k)) + F*(1+y)^(-k)
knitr::kable(tibble(`到期收益率（%）`= YTM, `半年收益率（%）`= YTM/m, `债券价格` = p), digits=2)| 到期收益率（%） | 半年收益率（%） | 债券价格 |
|---|---|---|
| 6 | 3.0 | 97.29 |
| 7 | 3.5 | 94.67 |
| 8 | 4.0 | 92.14 |
| 9 | 4.5 | 89.68 |
| 10 | 5.0 | 87.31 |
- 牛顿迭代法求解到期收益率
ytm.newton <- function(P=97.29, F=100, alpha=0.05, n=3, m=2){
k <- n*m
f <- function(y){
F*alpha/(m*y)*(1 - (1+y)^(-k)) + F*(1+y)^(-k) - P
}
x0 <- alpha/m + (F/P)^(1/k) - 1 #初始值，设定为：票面利率+零息债利率
eps <- 1E-6 # 方程左边绝对值小于此值时迭代结束
delta <- 1E-6 # 数值微分的步长
max.iter <- 100 #最大迭代次数
iter <- 0
repeat{
iter <- iter + 1
dfx <- (f(x0 + delta) - f(x0 - delta))/(2*delta) #用数值微分估计导数，可以根据泰勒公式推导
x0 <- x0 - f(x0)/dfx #牛顿迭代法的核心步骤
cat(iter, x0, f(x0), "\n")
if(iter >= max.iter || abs(f(x0))<eps) break
}
x0 # 结果是1/m年的到期收益率
}
ytm.newton(P=97.29, F=100, alpha=0.05, n=3, m=2)## 1 0.03000207 0.0003023458
## 2 0.03000264 5.709069e-10## [1] 0.03000264- 牛顿迭代法解方程的R代码
# 定义目标函数
f <- function(R) {
return(-100000*(1+R)^36*R+((1+R)^36-1)*3111.11)
}
# 定义牛顿迭代法函数
newton <- function(f, x0, eps, delta, max_iter) {
x <- x0
iter <- 0
while (iter < max_iter) {
# 计算函数值和导数值
fx <- f(x)
fpx <- (f(x+delta)-f(x-delta))/(2*delta)
# 计算新的近似解
x_new <- x - fx / fpx
# 判断是否满足停止条件
if (abs(f(x_new)) < eps) {
break
}
# 更新近似解和迭代次数
x <- x_new
iter <- iter + 1
cat(iter, x, f(x), "\n")
}
# 返回近似解和迭代次数
return(list(approximation = x, iterations = iter))
}
# 调用牛顿迭代法函数
result <- newton(f, x0 = 0.05, eps = 1e-12, delta = 1e-12, max_iter = 100)## 1 0.03527532 -4561.7
## 2 0.02383637 -1412.41
## 3 0.01570871 -411.6399
## 4 0.01053582 -110.5149
## 5 0.007697036 -25.39278
## 6 0.006518954 -3.783397
## 7 0.006270337 -0.159676
## 8 0.006258884 -0.0003612919
## 9 0.006258858 2.552359e-07
## 10 0.006258858 5.854872e-11ytm = result$approximation
# 打印结果
cat("月利率 =", ytm, "年利率 =", (1+ytm)^12-1)## 月利率 = 0.006258858 年利率 = 0.07774644- \(F=100\),\(P=95\), 到期持有天数在28到364之间变化时两种不同算法的变化曲线
F <- 100
P <- 95
curve((F-P)/F*360/x*100, 28, 364, col="red",
xlab="Days", ylab="YTM")
curve(((F/P)^(365.25/x)-1)*100, 28, 364,
add=TRUE, lty=2, col="blue")
legend("topright", lty=c(1,2), col=c("red", "blue"),
legend=c("Standard", "Exact"))
- 十年期国债
library(ggplot2)
file1 <- "D:/齐安静 教学/时间序列分析/美国十年期国债收益率历史数据.csv"
TNX <- read.table(file=file1, header=TRUE, sep=",")
colnames(TNX) <- c("date","close","open","high","low","increase")
TNX$date <- as.Date(TNX$date)
head(TNX)## date close open high low increase
## 1 2023-08-17 4.306 4.255 4.316 4.255 1.12%
## 2 2023-08-16 4.258 4.211 4.280 4.170 0.94%
## 3 2023-08-15 4.219 4.199 4.274 4.166 0.42%
## 4 2023-08-14 4.201 4.172 4.215 4.146 1.03%
## 5 2023-08-11 4.158 4.123 4.178 4.080 1.23%
## 6 2023-08-10 4.107 4.009 4.113 3.957 2.27%ggplot(data = TNX) +
geom_line(aes(x=date, y=close), color="green") +
labs(title ="十年期国债收盘价")+
theme_bw()
## 1.3 波动率 - 芝加哥期权交易所波动率指数
library(readr)
library(lubridate)
library(xts)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric##
## ######################### 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: 'xts'## The following objects are masked from 'package:dplyr':
##
## first, lastlibrary(quantmod)## Loading required package: TTR## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zootmp.s <- read_file("D:/齐安静 教学/时间序列分析/北大/ftsdata/d-vix0411.txt")
tmp.s <- gsub("\t", " ", tmp.s, fixed=TRUE)
d.vix <- read_table(tmp.s,
col_types=cols(.default = col_double()))
rm(tmp.s)
vix <- xts(d.vix[4:7],
make_date(d.vix[["year"]], d.vix[["mon"]], d.vix[["day"]]))
chartSeries(
vix, type="line", TA=NULL,
major.ticks="years", minor.ticks=FALSE,
theme="white", name="CBOE Volatility Index"
)
## 1.4 收益率的分布性质 - 苹果公司日数据
library(dplyr)
aapl <- read.csv("D:/齐安静 教学/时间序列分析/aapl.csv")
aapl$date <- as.Date(aapl$Date)
aapl <- mutate(aapl, logret = log(aapl$Adj.Close))
aapl <- mutate(aapl, diff1 = aapl$logret - lag(aapl$logret))
head(aapl)## Date Open High Low Close Adj.Close Volume date
## 1 2001-08-17 0.321429 0.329464 0.321250 0.322679 0.273889 208426400 2001-08-17
## 2 2001-08-20 0.323929 0.325536 0.318036 0.323571 0.274646 252302400 2001-08-20
## 3 2001-08-21 0.323929 0.323929 0.316071 0.320000 0.271615 185701600 2001-08-21
## 4 2001-08-22 0.320357 0.325893 0.314464 0.325179 0.276011 173975200 2001-08-22
## 5 2001-08-23 0.325000 0.327500 0.313929 0.318036 0.269948 217078400 2001-08-23
## 6 2001-08-24 0.321429 0.332500 0.315179 0.331607 0.281467 290332000 2001-08-24
## logret diff1
## 1 -1.295032 NA
## 2 -1.292272 0.002760081
## 3 -1.303370 -0.011097373
## 4 -1.287315 0.016055098
## 5 -1.309526 -0.022211372
## 6 -1.267740 0.041785864tail(aapl)## Date Open High Low Close Adj.Close Volume date
## 5528 2023-08-09 180.87 180.93 177.01 178.19 177.9497 60378500 2023-08-09
## 5529 2023-08-10 179.48 180.75 177.60 177.97 177.7300 54686900 2023-08-10
## 5530 2023-08-11 177.32 178.62 176.55 177.79 177.7900 51988100 2023-08-11
## 5531 2023-08-14 177.97 179.69 177.31 179.46 179.4600 43675600 2023-08-14
## 5532 2023-08-15 178.88 179.48 177.05 177.45 177.4500 43622600 2023-08-15
## 5533 2023-08-16 177.13 178.54 176.50 176.57 176.5700 46876800 2023-08-16
## logret diff1
## 5528 5.181501 -0.0089946703
## 5529 5.180266 -0.0012354432
## 5530 5.180603 0.0003375169
## 5531 5.189952 0.0093493417
## 5532 5.178689 -0.0112635187
## 5533 5.173717 -0.0049714242ggplot(data = aapl) +
geom_line(aes(x=date, y=diff1), color="green") +
labs(title ="苹果公司对数收益率")+
theme_bw()## Warning: Removed 1 row containing missing values (`geom_line()`).
- 2017年以来的数据
newdata <- aapl %>%
filter(aapl$Date >= "2017-01-01")
# newdata <- subset(aapl, date >= 2017-01-01)
head(newdata)## Date Open High Low Close Adj.Close Volume date
## 1 2017-01-03 28.9500 29.0825 28.6900 29.0375 27.05931 115127600 2017-01-03
## 2 2017-01-04 28.9625 29.1275 28.9375 29.0050 27.02902 84472400 2017-01-04
## 3 2017-01-05 28.9800 29.2150 28.9525 29.1525 27.16647 88774400 2017-01-05
## 4 2017-01-06 29.1950 29.5400 29.1175 29.4775 27.46934 127007600 2017-01-06
## 5 2017-01-09 29.4875 29.8575 29.4850 29.7475 27.72094 134247600 2017-01-09
## 6 2017-01-10 29.6925 29.8450 29.5750 29.7775 27.74890 97848400 2017-01-10
## logret diff1
## 1 3.298031 0.002845062
## 2 3.296911 -0.001119983
## 3 3.301984 0.005072499
## 4 3.313070 0.011086725
## 5 3.322188 0.009117684
## 6 3.323196 0.001008043ggplot(data = newdata) +
geom_line(aes(x=date, y=diff1), color="green") +
labs(title ="2017年以来苹果公司对数收益率")+
theme_bw()
- 直方图
x <- coredata(aapl$diff1[-1])
mean1=mean(x,na.rm=TRUE)
sd1=sd(x,na.rm=TRUE)
hist(x,xlab="对数收益率",ylab = "频率",freq = FALSE,cex.lab=0.7,ylim = c(0,28),xlim=c(-0.15,0.15),main="对数收益率直方图")
d <- seq(from=min(x),to=max(x),by=0.001)
lines(x=d,y=dnorm(d,mean1,sd1),lty=2,col=2)
lines(density(x),lty=4)
- 正态QQ图
x <- coredata(aapl$diff1)
qqnorm(x, main="Apple Log Returns")
qqline(x, col="red")
- 十年期美国国债
- 对数收益率
library(ggplot2)
file1 <- "D:/齐安静 教学/时间序列分析/^TNX.csv"
TNX <- read.table(file=file1, header=TRUE, sep=",", na="null")
rows <- which(is.na(TNX$Close))
TNX <- TNX[-rows,]
TNX$date <- as.Date(TNX$Date)
TNX <- mutate(TNX, logret = log(TNX$Adj.Close))
TNX <- mutate(TNX, diff1 = TNX$logret - lag(TNX$logret))
head(TNX)## Date Open High Low Close Adj.Close Volume date logret
## 1 2018-08-20 2.853 2.853 2.821 2.823 2.823 0 2018-08-20 1.037800
## 2 2018-08-21 2.844 2.853 2.833 2.844 2.844 0 2018-08-21 1.045212
## 3 2018-08-22 2.819 2.835 2.808 2.823 2.823 0 2018-08-22 1.037800
## 4 2018-08-23 2.821 2.826 2.814 2.821 2.821 0 2018-08-23 1.037091
## 5 2018-08-24 2.832 2.850 2.817 2.826 2.826 0 2018-08-24 1.038862
## 7 2018-08-27 2.830 2.851 2.828 2.848 2.848 0 2018-08-27 1.046617
## diff1
## 1 NA
## 2 0.0074113627
## 3 -0.0074113627
## 4 -0.0007087173
## 5 0.0017708522
## 7 0.0077547093ggplot(data = TNX) +
geom_line(aes(x=date, y=diff1), color="green") +
labs(title ="US 10 Years T Notes Log Return")+
theme_bw()## Warning: Removed 1 row containing missing values (`geom_line()`).
-直方图
x <- coredata(TNX$diff1[-1])
mean1=mean(x,na.rm=TRUE)
sd1=sd(x,na.rm=TRUE)
hist(x,xlab="对数收益率",ylab = "频率",freq = FALSE,cex.lab=0.7,ylim = c(0,16),xlim=c(-0.15,0.15),main="TNX对数收益率直方图")
d <- seq(from=min(x),to=max(x),by=0.001)
lines(x=d,y=dnorm(d,mean1,sd1),lty=2,col=2)
lines(density(x),lty=4)
- 正态QQ图
x <- coredata(TNX$diff1)
qqnorm(x, main="Q-Q Plot of US 10 Years T-Notes Log Returns")
qqline(x, col="red")
- 欧元对美元的日汇率
library(quantmod)
DEXUSEU <- getSymbols("DEXUSEU", src="FRED",
auto.assign=FALSE)logret.DEXUSEU <- diff(log(DEXUSEU))*100 # 对数收益率chartSeries(
DEXUSEU, type="l", TA=NULL,
subset="2012/2022",
name="Euro to USD Price",
theme="white", major.ticks="years",
minor.ticks=FALSE)
- 对数收益率
chartSeries(
logret.DEXUSEU, type="l", TA=NULL,
subset="2012/2022",
name="Euro to USD Log Returns",
theme="white", major.ticks="years", minor.ticks=FALSE)
- 直方图
x <- coredata(logret.DEXUSEU)
mean1=mean(x,na.rm=TRUE)
sd1=sd(x,na.rm=TRUE)
rows <- which(is.na(x))
x <- x[-rows,]
hist(x, main="Euro to USD Log Returns", xlab="", ylab="",freq=FALSE,ylim=c(0,0.9),xlim=c(-2,2))
print(min(x))## [1] -3.003101d <- seq(from=min(x),to=max(x),by=0.001)
lines(x=d,y=dnorm(d,mean1,sd1),lty=2,col=2)
lines(density(x),lty=4)
- Q-Q图
x <- coredata(logret.DEXUSEU)
qqnorm(x, main="Euro to USD Log Returns")
qqline(x, col="red")
1.5.2 收益率分布案例
- 苹果公司2011-2017对数收益率
library(quantmod)
AAPL <- getSymbols("AAPL", src="tiingo", auto.assign=FALSE, api.key="9667b4325f683d9d0ae0734af175d886eceb7781")library(fBasics)## Warning: package 'fBasics' was built under R version 4.3.1##
## Attaching package: 'fBasics'## The following object is masked from 'package:TTR':
##
## volatilityx <- coredata(diff(
log(AAPL["2011/2017"]$AAPL.Adjusted)))[-1]
fBasics::basicStats(x)## x
## nobs 1760.000000
## NAs 0.000000
## Minimum -0.131875
## Maximum 0.085022
## 1. Quartile -0.006914
## 3. Quartile 0.009382
## Mean 0.000789
## Median 0.000574
## Sum 1.388788
## SE Mean 0.000377
## LCL Mean 0.000050
## UCL Mean 0.001528
## Variance 0.000250
## Stdev 0.015809
## Skewness -0.323132
## Kurtosis 5.490065- 直方图
hist(x, nclass=30, main="Apple Log Returns")
- 作核密度估计图并叠加正态密度估计
tmp.1 <- density(x, na.rm=TRUE)
tmp.x <- seq(min(x, na.rm=TRUE), max(x, na.rm=TRUE),
length.out=100)
tmp.y <- dnorm(tmp.x, mean(x, na.rm=TRUE),
sd(x, na.rm=TRUE))
tmp.ra <- range(c(tmp.1$y, tmp.y), na.rm=TRUE)
plot(tmp.1, main="Apple Log Return",
ylim=tmp.ra)
lines(tmp.x, tmp.y, lwd=2, col="red")
legend("topleft", lwd=c(1,2),
col=c("black", "red"),
legend=c("Kernel density Est.",
"Parametric normal density est."))
- t检验 `
t.test(x)##
## One Sample t-test
##
## data: x
## t = 2.094, df = 1759, p-value = 0.03641
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 4.999128e-05 1.528177e-03
## sample estimates:
## mean of x
## 0.000789084- 偏度
tmp <- fBasics::basicStats(x)["Skewness", 1]; tmp## [1] -0.323132tmp/sqrt(6/length(x))## [1] -5.534274- 峰度
tmp <- fBasics::basicStats(x)["Kurtosis", 1]; tmp## [1] 5.490065tmp/sqrt(24/length(x))## [1] 47.01411- JB检验
fBasics::normalTest(x, method="jb")##
## Title:
## Jarque - Bera Normalality Test
##
## Test Results:
## STATISTIC:
## X-squared: 2248.7883
## P VALUE:
## Asymptotic p Value: < 2.2e-161.6 可视化
- 时间序列曲线
chartSeries(
AAPL, type="line",
subset="2017", TA=NULL,
theme="white", name="Apple Stock Price",
major.ticks="months", minor.ticks=FALSE)
- 带有7日均线
chartSeries(
AAPL, type="line",
subset="2017", TA="addSMA(7)",
theme="white", name="Apple Stock Price",
major.ticks="months", minor.ticks=FALSE)
- 用ggplot2作图
library(ggplot2)
p <- ggplot(data = tibble(
time=index(AAPL),
price=coredata(AAPL)[,"AAPL.Adjusted",drop=TRUE]),
mapping = aes(x=time, y=price))
p + geom_line()
- 装饰
p + geom_line() +
scale_x_date(
name = "year",
date_breaks = "5 years",
date_labels = "%Y", #表示用年来做刻度
date_minor_breaks = "1 year",
expand = c(0, 0)) 
- 人为指定年
p + geom_line() +
scale_x_date(
name = "year",
breaks = lubridate::make_date(c(2008, 2011, 2014, 2017, 2020,2023)),
date_labels = "%Y",
date_minor_breaks = "1 year",
expand = c(0, 0)) 
- 日期旋转
p + geom_line() +
scale_x_date(
name = "year",
date_labels = "%Y-%m-%d",
expand = c(0, 0)) +
theme(
axis.text.x = element_text(
angle = 45, vjust = 1, hjust = 1) )
### 1.6.2 k线图
- AAPL 2017k线图
chartSeries(
AAPL,
subset="2017",
theme="white", name="Apple Stock Price",
major.ticks="months", minor.ticks=FALSE)
- 仅画最后一个月
chartSeries(
AAPL,
subset="2017-12",
theme="white", name="Apple Stock Price",
major.ticks="auto", minor.ticks=FALSE)
### 1.6.3 密度估计 - 苹果公司2017年简单收益率的直方图
x <- simple.return(AAPL["2017","AAPL.Adjusted"])
hist(x, main="Apple Stock 2017 Simple Return", xlab="")
- 核密度线
tmp.1 <- density(x, na.rm=TRUE)
tmp.x <- seq(min(x, na.rm=TRUE), max(x, na.rm=TRUE),
length.out=100)
tmp.y <- dnorm(tmp.x, mean(x, na.rm=TRUE),
sd(x, na.rm=TRUE))
tmp.ra <- range(c(tmp.1$y, tmp.y), na.rm=TRUE)
plot(tmp.1, main="Apple Stock 2017 Simple Return",
ylim=tmp.ra)
lines(tmp.x, tmp.y, lwd=2, col="red")
legend("topright", lwd=c(1,2),
col=c("black", "red"),
legend=c("Kernel density Est.",
"Parametric normal density est."))
### 1.6.4 q-q图
qqnorm(x, main="Apple Stock 2017 Simple Return")
qqline(x, col="red")
### 1.6.5 散点图与回归直线 -
读入IBM和S&P的月度收益率数据，从1926-01到2011-09,
并转换为xts格式，时间下标改为yearmon格式：
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.0432ibmsp <- xts(as.matrix(d[,2:3]), d$date)
tclass(ibmsp) <- "yearmon"
head(ibmsp)## ibm sp
## Jan 1926 -0.010381 0.022472
## Feb 1926 -0.024476 -0.043956
## Mar 1926 -0.115591 -0.059113
## Apr 1926 0.089783 0.022688
## May 1926 0.036932 0.007679
## Jun 1926 0.068493 0.043184- 序列图
chartSeries(
ibmsp[,"ibm"], type="lines",
theme="white", name="IBM returns",
major.ticks="years", minor.ticsk=FALSE)
- S&P标普指数
chartSeries(
ibmsp[,"sp"], type="lines",
theme="white", name="S&P returns",
major.ticks="years", minor.ticsk=FALSE)
- 散点图
d <- as.matrix(coredata(ibmsp))
plot(d[,"sp"], d[,"ibm"], pch=16, cex=0.2,
xlab="S&P", ylab="IBM")
- 相关系数
d <- as.matrix(coredata(ibmsp))
cor.test(d[,"sp"], d[,"ibm"])##
## Pearson's product-moment correlation
##
## data: d[, "sp"] and d[, "ibm"]
## t = 26.664, df = 1027, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6020164 0.6743519
## sample estimates:
## cor
## 0.6395979- 散点图+回归线
d <- as.data.frame(coredata(ibmsp))
plot(d[["sp"]], d[["ibm"]], pch=16, cex=0.2,
xlab="S&P", ylab="IBM")
lm1 <- lm(ibm ~ sp, data=d)
print(summary(lm1))##
## Call:
## lm(formula = ibm ~ sp, data = d)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.279155 -0.032137 -0.002261 0.030647 0.280313
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.008974 0.001706 5.259 1.76e-07 ***
## sp 0.818776 0.030707 26.664 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05443 on 1027 degrees of freedom
## Multiple R-squared: 0.4091, Adjusted R-squared: 0.4085
## F-statistic: 711 on 1 and 1027 DF, p-value: < 2.2e-16abline(lm1, col="red", lwd=2)