Homework 2

学号：1940802 姓名：刘睿奇

首先导入基本的模块。

library(fBasics)
library(TTR)
library(xts)
library(zoo)
library(quantmod)
library(moments)

1、 下表为14名大学生的英语四级成绩和六级考试成绩。试求：

序号 1 2 3 4 5 6 7 8 9 10 11 12 13 14 四级考 98 85 89 84 81 70 92 67 84 80 60 81 73 70 六级考 90 82 88 80 82 66 88 68 84 77 64 79 75 73

先在Excel中将原表的行列转置如下（并改中文为英文方便识别）： Num CET4 CET6 1 98 90 2 85 82 3 89 88 4 84 80 5 81 82 6 70 66 7 92 88 8 67 68 9 84 84 10 80 77 11 60 64 12 81 79 13 73 75 14 70 73

读取数据：

da = read.table("clipboard",header = T)  #需要先复制上面的表格
print(da)

a）求出这两者之间的相关系数，画出它们的散点图；

cor(da$CET4,da$CET6)  #求相关系数

[1] 0.9641479

plot(da$CET4, da$CET6)

b）求出各考试的均值、标准差、偏度、峰度；

basicStats(da[,2:3])

以上的Mean为均值，Stdev为标准差，Skewness为偏度，Kurtosis为（超额）峰度。

c）画出各考试的盒形图（在一张图上的对比盒形图）；

boxplot(da[,2:3])
title(main = 'CET4 (left) & CET6 (right)')

d）画出各考试的茎叶图。

print('CET4 STEM')

[1] "CET4 STEM"

stem(da$CET4)

  The decimal point is 1 digit(s) to the right of the |

  6 | 07
  7 | 003
  8 | 0114459
  9 | 28

print('CET6 STEM')

[1] "CET6 STEM"

stem(da$CET6)

  The decimal point is 1 digit(s) to the right of the |

  6 | 468
  7 | 3579
  8 | 022488
  9 | 0

2、某种电子元件的平均寿命x(单位：小时)服从正态分布，现测得16只元件的寿命分别为159、280、101、212、224、379、179、264、222、362、168、149、260、485、170，试求元件寿命的极大似然估计及区间估计，并问有否理由认为元件的平均寿命地大于225小时（α=0.05）。

x <- c(159,280,101,212,224,379,179,264,222,362,168,149,260,485,170)
hist(x, breaks = 50*(0:12)，probability = TRUE)
d1 <- density(x)
lines(d1)  #做密度曲线

1)极大似然估计 定义似然函数

LL <- function(params, data) {
  R <-  dnorm(data, params[1], params[2])
  ll <- sum(log(R))
  return(-ll)
}

参数估计

da.res1 <- nlminb(c(293, 50),
                 LL,
                 data = x,
                 lower = c(0.0001, 0.0001),
                 upper = c(Inf, Inf)
                 )
da.res1$par  #估计结果

[1] 240.93333  98.69986

由上可知，极大似然估计的结果是：mean == 240.93333，sd == 98.69986

做以上估计结果的拟合图如下：

xl <- seq(0,600,1)
f <- dnorm(xl,da.res1$par[1],da.res1$par[2])
hist(x, breaks = 50*(0:12)，probability = TRUE, col = 'lightblue',)
lines(xl, f, col = "darkblue", lwd = 2)

lines(d1, col = 'red', lwd = 2)

legend('topright',lwd = 2, legend = c('res','density'), col = c('darkblue','red')) 

2）区间估计

t.test(x, conf.level = 0.95) 

    One Sample t-test

data:  x
t = 9.1337, df = 14, p-value = 2.835e-07
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 184.3568 297.5099
sample estimates:
mean of x 
 240.9333 

由上面的结果可知，95%的置信区间为184.3568到297.5099。

3）是否有理由认为元件的平均寿命大于225小时（α=0.05）

H0: 平均寿命等于225 Ha: 平均寿命大于225 在t.test函数中，将alternative设置为大于（greater），mu设置为225，进行检验

t.test(x, alternative = 'greater', mu = 225, conf.level = 0.95)

    One Sample t-test

data:  x
t = 0.60402, df = 14, p-value = 0.2777
alternative hypothesis: true mean is greater than 225
95 percent confidence interval:
 194.4724      Inf
sample estimates:
mean of x 
 240.9333 

由上可知，得到的P值为0.2777（>0.05），显然无法拒绝原假设，即没有理由认为原件的平均寿命大于225小时。

3. 考虑从2001年9月1日到2011年9月30日美国运通公司（AXP）、CRSP价值权重指数（VW）、CRSO的等权重指数（EW）以及S&P综合指数的日简单收益率。收益率中包含有支付的股息。数据来自d-axp3dx-0111.txt.对AXP股票的对数收益的均值等于零的原假设进行检验。在5%的显著性水平下，得出你的结论。

读取数据：

da <- read.table('d-axp3dx-0111.txt', header = T, row.names = 1)
head(da)

daAXP <- da[1]
head(daAXP)

计算对数收益率，并进行t检验

daAXP_L <- log(daAXP + 1)
head(daAXP_L)

t.test(daAXP_L$axp)

    One Sample t-test

data:  daAXP_L$axp
t = 0.35999, df = 2534, p-value = 0.7189
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -0.0008360686  0.0012120714
sample estimates:
   mean of x 
0.0001880014 

P值的结果为0.7189（>0.05），故无法拒绝原假设，AXP股票的对数收益的均值等于零。

4. 考虑从1940年1月到2011年9月S&P综合指数的月股票收益率（m-ge3dx-4011.txt），进行下面的检验，在5%的显著性水平下得出你的结论。

读取数据：

da <- read.table('m-ge3dx-4011.txt', header = T, row.names = 1)
head(da)

可见sp是fctr属性，于是按照老师的说明，进行数据处理：

sp.num <- as.numeric(as.character(da$sp))

NAs introduced by coercion

da$sp <- sp.num

#填充缺失值为非缺失值的平均数
da$sp[is.na(da$sp)] <- mean(da$sp[!is.na(da$sp)])

head(da)

可见sp已经转换为dbl属性了。

1. 检验假设H0: μ=0，其备择假设为H1：μ<>0，这里μ是收益率的均值；

t.test(da$sp)

    One Sample t-test

data:  da$sp
t = 4.2488, df = 860, p-value = 2.384e-05
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 0.003323913 0.009031515
sample estimates:
  mean of x 
0.006177714 

P值很小，可以拒绝原假设，可知μ不等于0。

2. 检验假设H0: m3=0，其备择假设为H1：m3<>0，这里m3是收益率的偏度；

m3 <- skewness(da$sp)  #计算偏度
m3

[1] -0.5908271

agostino.test(da$sp)

    D'Agostino skewness test

data:  da$sp
skew = -0.59083, z = -6.63877, p-value = 3.163e-11
alternative hypothesis: data have a skewness

从检验的结果可知P值很小，可以拒绝原假设，可知偏度不等于0。

3. 检验假设H0: K=3，其备择假设为H1：K<>0，这里K是收益率的峰度；

k <- kurtosis(da$sp)+3  #计算峰度
k

[1] 8.380011

anscombe.test(da$sp) 

    Anscombe-Glynn kurtosis test

data:  da$sp
kurt = 5.3800, z = 7.2265, p-value = 4.955e-13
alternative hypothesis: kurtosis is not equal to 3

从检验的结果可知P值很小，可以拒绝原假设，可知峰度不等于3。

5. 从芝加哥的联邦储备银行得到日汇率，数据是经过纽约联邦储备银行认证的纽约市每日中午买入价。考虑从2007年1月2日到2011年11月30日美元对英镑、日元的汇率（数据见：d-fx-usjp-0711.txt）。 计算：

读取数据：

da <- read.table('d-fx-usjp-0711.txt', header = T)
length(da$rate)

[1] 1238

head(da)

tail(da)

1. 每个汇率的日对数收益率；

rate.lrr <- diff(log(da$rate))
length(rate.lrr)

[1] 1237

head(rate.lrr)

[1]  0.0062917032 -0.0025119330 -0.0042848207 -0.0023602810  0.0080693049
[6]  0.0005858476

tdx <- da$year + da$mon/12 + da$day/12/31  #产生简易的时间指数
plot(tdx[2:nrow(da)], rate.lrr,type = 'l',xlab = 'year')

2. 汇率的日对数收益率的样本均值、标准差、偏度、超额峰度、最大值和最小值；

basicStats(rate.lrr)

以上数据中，Mean为样本均值，Stdev为标准差，Skewness为偏度，Kurtosis为超额峰度，Maximum为最大值，Minimum为最小值。

3. 将作图区间分为两个，在其中画出美元/日元汇率的走势图及日对数收益率的密度函数（包括核密度估计和正态密度估计）；

par(mfrow = c(1,2))
#par(mfcol = c(2,1))

plot(tdx, da$rate,type = 'l',xlab = 'year')
title(main = '走势图')

d2 <- density(rate.lrr)
plot(d2, col = 'red', lwd = 2,xlim = c(-0.06,0.06), ylim = c(0,80),title ='')
x <- seq(-0.1,0.1,0.001)
y <- dnorm(x, mean(rate.lrr), stdev(rate.lrr))
lines(x, y, col = "darkblue", lwd = 2)

legend('topright',lwd = 2, legend = c('norm','dens'), col = c('darkblue','red')) 

4. 检验假设H0: μ=0，其备择假设为H1：μ<>0，这里μ表示美元/日元汇率的日对数收益率的均值，在5%的显著性水平下给出你的结论。

t.test(rate.lrr)

    One Sample t-test

data:  rate.lrr
t = -1.6133, df = 1236, p-value = 0.1069
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -7.638670e-04  7.448262e-05
sample estimates:
    mean of x 
-0.0003446922 

p-value == 0.1069 > 0.05，故无法拒绝原假设，可得μ == 0

6. 请下载苹果公司（AAPL）的股票价格（从2019年初至今），试求：

读取数据：

getSymbols('AAPL', from = '2019-01-01')

[1] "AAPL"

length(AAPL)

[1] 2298

head(AAPL)

           AAPL.Open AAPL.High AAPL.Low AAPL.Close AAPL.Volume AAPL.Adjusted
2019-01-02    154.89    158.85   154.23     157.92    37039700      154.7950
2019-01-03    143.98    145.72   142.00     142.19    91312200      139.3763
2019-01-04    144.53    148.55   143.80     148.26    58607100      145.3261
2019-01-07    148.70    148.83   145.90     147.93    54777800      145.0027
2019-01-08    149.56    151.82   148.52     150.75    41025300      147.7669
2019-01-09    151.29    154.53   149.63     153.31    45099100      150.2762

tail(AAPL)

           AAPL.Open AAPL.High AAPL.Low AAPL.Close AAPL.Volume AAPL.Adjusted
2020-07-01    365.12    367.36   363.91     364.11    27684300        364.11
2020-07-02    367.85    370.47   363.64     364.11    28510400        364.11
2020-07-06    370.00    375.78   369.87     373.85    29663900        373.85
2020-07-07    375.41    378.62   372.23     372.69    28106100        372.69
2020-07-08    376.72    381.50   376.36     381.37    29273000        381.37
2020-07-09    385.05    385.27   378.69     383.01    31362900        383.01

1. 画出其收盘价的时间序列图及五日均线、十日均线、三十日均线；

source('ma.R')
x1 <- as.numeric(AAPL$AAPL.Close)
ma(x1,5)

ma(x1,10)

ma(x1,30)

2. 试对AAPL股票价画其烛线图（股票条形图）；

X <- AAPL[,1:4]
x2 <- cbind(as.numeric(X[,1]),as.numeric(X[,2]),as.numeric(X[,3]),as.numeric(X[,4]))

source('ohlc.R')
ohlc_plot(x2, xl = 'days', yl = 'price', title = 'Apple Stock')

3. 计算收盘价的对数收益率，并画出其收益率图（c与d的图一个在上，一个在下）；
4. 求出其ACF图（lag=24）。

par(mfcol = c(2,1))

AAPL.lrr <- diff(log(AAPL$AAPL.Close))
plot(AAPL.lrr)

acf(AAPL.lrr$AAPL.Close[2:nrow(AAPL.lrr)], lag = 24)

上面的ACF稍显不清晰，做大图如下：

acf(AAPL.lrr$AAPL.Close[2:nrow(AAPL.lrr)], lag = 24)