Homework 3

学号：1940802 姓名：刘睿奇

首先导入基本的模块。

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

Attaching package: 㤼㸱TSA㤼㸲

The following objects are masked from 㤼㸱package:moments㤼㸲:

    kurtosis, skewness

The following objects are masked from 㤼㸱package:timeDate㤼㸲:

    kurtosis, skewness

The following objects are masked from 㤼㸱package:stats㤼㸲:

    acf, arima

The following object is masked from 㤼㸱package:utils㤼㸲:

    tar

1. 从正态分布N(100,100)中随机产生1000个随机数,

1. 作出这1000个正态随机数的直方图;

da <- rnorm(1000, 100, sqrt(100))
hist(da,probability = T)
lines(density(da))

2. 从这1000个随机数中随机有放回地抽取500个, 作出其直方图;

da500 <- sample(da, size = 500, replace = T)
hist(da500,probability = T)
lines(density(da500))

3. 比较它们的样本均值与样本方差.

print('Mean')

[1] "Mean"

mean(da)

[1] 100.2

mean(da500)

[1] 99.75

print('Var')

[1] "Var"

var(da)

[1] 108.3

var(da500)

[1] 105.2

2.假定某校100名女生的血清总蛋白含量(g/L)服从均值为75, 标准差为9,并假定数据由下面的命令产生

options(digits=4)
da <- rnorm(100,75,9)

根据产生的数据: (a) 计算样本均值、方差、标准差、极差、四分位极差、变异系数、偏度、峰度和五数概括;

bda <- basicStats(da)
bda

如上所示，Mean为样本均值，Variance为方差，Stdev为标准差，Skewness为偏度，Kurtosis为（超额）峰度

dra <- bda["Maximum",1] - bda["Minimum",1]
print('dra')

[1] "dra"

dra

[1] 36.02

print('IQR')

[1] "IQR"

IQR(da)

[1] 10.52

print('cv')

[1] "cv"

cv <- bda["Mean",1]/bda["Stdev",1]

dra求得的结果为极差，IQR求得的结果为四分位极差，cv求得的结果为变异系数

summary(da)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   57.3    70.8    76.0    76.0    81.3    93.3 

以上结果除了Mean剩下的五个数即为五数概括

2. 画出直方图、核密度估计曲线、QQ图;

hist(da,probability = T)
lines(density(da))

qqnorm(da, main="Normality Check via QQ Plot")   #QQ图
qqline(da, col='red')   

3. 画出茎叶图、盒形图

stem(da)

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

  5 | 7
  6 | 0022344
  6 | 5666788889999
  7 | 000011223334444444
  7 | 5555555566667777888888889999
  8 | 0000011112223333344444
  8 | 5677777799
  9 | 3

boxplot(da)
title(main = 'Boxplot')

3. 考虑美国从1947年第1季度到2011年第3季度的季度实际GNP，该数据存放于文件q-GNPC96.txt中，数据已做了季节调整，以2005年GDP为基础进行了通胀调整，以10亿美元为单位。假设x_t代表GNP增长率的时间序列数据.

读取数据：

da <- read.table('q-GNPC96.txt', header = T)
length(da$gnp)

[1] 259

head(da)

tail(da)

计算增长率x_t（简单收益率）：

x_t <- 1 / (da$gnp[2:nrow(da)] / diff(da$gnp[2:nrow(da)]) - 1)

longer object length is not a multiple of shorter object length

length(x_t)

[1] 258

head(x_t)

[1] -0.0008429  0.0156643  0.0165671  0.0190084  0.0053833  0.0015469

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

1. 通过ar命令，应用AIC准则，可以为x_t识别一个几阶AR()模型？拟合这个模型，拟合的模型充分吗？为什么？

模型定阶：

mm11 <- ar(x_t,method='mle') #用AIC准则，自动为AR(P)定阶，方法为极大似然估计
mm11$order #寻找合适的阶p

[1] 3

print(mm11$aic) #查看mm1中aic的值

      0       1       2       3       4       5       6       7       8       9 
38.5705  3.1667  2.4577  0.0000  0.1469  0.3142  1.9640  3.8857  5.8843  5.7963 
     10      11      12 
 7.7942  9.7707  4.9022 

aic <- mm11$aic
plot(c(0:(length(aic)-1)),aic,type='h',xlab='order',ylab='aic') #画aic竖线图
lines(0:(length(aic)-1),aic,lty=2) #画aic连线图

由上述结果可知，3的AIC最小（为0.000），可以识别为一个AR(3)的模型。

参数估计：

mm12 <- arima(x_t,order=c(3,0,0))
mm12

Call:
arima(x = x_t, order = c(3, 0, 0))

Coefficients:
        ar1    ar2     ar3  intercept
      0.344  0.145  -0.132      0.008
s.e.  0.062  0.065   0.062      0.001

sigma^2 estimated as 8.81e-05:  log likelihood = 838.2,  aic = -1668

模型检验： H0：残差为白噪声 Ha：残差不为白噪声

Box.test(mm12$residuals,lag=12,type='Ljung') 

    Box-Ljung test

data:  mm12$residuals
X-squared = 11, df = 12, p-value = 0.5

pv=1-pchisq(11,9) #去掉3个参数后，用自由度9检验p值
pv

[1] 0.2757

p-value > 0.05，无法拒绝原假设，故残差为白噪声，故模型是充分的。

2. 数据x_t的样本PACF识别的是几阶AR()时间序列模型？拟合这个模型，拟合的模型充分吗？ 为什么？

模型定阶：

pacf(x_t)

可以看到，PACF在3和12时还是显著的，但4~11都不显著，所以还是定为AR(3)模型更为合适。 参数估计、模型检验同第一问，不再重新计算。

4. 以 NYSE/AMEX/NASDAQ 的市场资本为基础考虑 CRSP Decile 1、2、5、9、10投资组合的月简单收益率. 该数据的时间区间是从1961年1月到2011年9月.

读取数据：

da <- read.table('m-dec125910-6111.txt', header=T, row.names = 1)
nrow(da)

[1] 609

head(da)

tail(da)

1. 对于 Decile 2 和 Decile 10的收益序列，在5%的显著性水平下检验：原假设是滞后阶数为1-12的自相关系数均为0. 给出你的结论.

par(mfcol=c(2,1))
acf(da$dec2,lag = 12)
acf(da$dec10,lag = 12)

对于Decile 2： H0：滞后阶数为1-12的自相关系数均为0 Ha：滞后阶数为1-12的自相关系数不全为0

Box.test(da$dec2,lag=12,type='Ljung')

    Box-Ljung test

data:  da$dec2
X-squared = 17, df = 12, p-value = 0.1

p-value > 0.05，故无法拒绝原假设。所以对于Decile 2，滞后阶数为1-12的自相关系数均为0

对于Decile 10： H0：滞后阶数为1-12的自相关系数均为0 Ha：滞后阶数为1-12的自相关系数不全为0

Box.test(da$dec10,lag=12,type='Ljung')

    Box-Ljung test

data:  da$dec10
X-squared = 48, df = 12, p-value = 4e-06

p-value < 0.05，故拒绝原假设。所以对于Decile 2，滞后阶数为1-12的自相关系数不全为0

2. 对于 Decile 2 的收益率序列建立一个 ARMA 模型, 对模型进行检验并写出拟合的模型.

模型定阶：

mm <- eacf(da$dec2,9,9)

AR/MA
  0 1 2 3 4 5 6 7 8 9
0 x o o o o o o o o o
1 x o o o o o o o o o
2 x x o o o o o o o o
3 x o x o o o o o o o
4 x o x o o o o o o o
5 x x x x x o o o o o
6 x x o x o x o o o o
7 x x o x x x x o o o
8 x x x x x o o x o o
9 x o o o x o x x x o

print(mm$eacf,digits=2)

       [,1]   [,2]     [,3]    [,4]    [,5]   [,6]    [,7]    [,8]    [,9]   [,10]
 [1,]  0.11 -0.054  0.00042  0.0332  0.0354 -0.064 -0.0088 -0.0627 -0.0352  0.0070
 [2,]  0.43 -0.053 -0.00083  0.0170  0.0211 -0.071  0.0111 -0.0476 -0.0431 -0.0115
 [3,]  0.21  0.194  0.06062 -0.0129  0.0172 -0.051 -0.0058 -0.0629 -0.0464 -0.0103
 [4,] -0.42  0.065  0.15791  0.0238 -0.0053 -0.055  0.0088 -0.0814 -0.0474 -0.0087
 [5,] -0.50  0.075  0.26943  0.0537 -0.0116 -0.069 -0.0078 -0.0415 -0.0367  0.0094
 [6,]  0.36  0.368  0.30068 -0.1653  0.2682 -0.074  0.0189 -0.0014 -0.0217 -0.0347
 [7,]  0.16 -0.472 -0.03422  0.3843  0.0681 -0.215  0.0292 -0.0095 -0.0034 -0.0428
 [8,]  0.14 -0.459 -0.04374  0.3828  0.1496 -0.110 -0.1439 -0.0247  0.0064 -0.0408
 [9,] -0.26  0.160 -0.13143  0.2240 -0.2530  0.026  0.0390 -0.1192 -0.0546 -0.0425
[10,]  0.32 -0.073  0.04478 -0.0072 -0.2709  0.019  0.1002 -0.1569  0.1428  0.0346

可见三角形的顶点在（0，0），虽然这个位置不为“o”，但是三角形的顶点应该在这里。选用ARMA(0,0)模型。

参数估计：

m1 <- arima(da$dec2, order = c(0,0,0))
m1

Call:
arima(x = da$dec2, order = c(0, 0, 0))

Coefficients:
      intercept
          0.009
s.e.      0.002

sigma^2 estimated as 0.00226:  log likelihood = 991.4,  aic = -1981

模型检验：

tsdiag(m1) 

Box.test(m1$residuals,lag=12,type='Ljung') 

    Box-Ljung test

data:  m1$residuals
X-squared = 17, df = 12, p-value = 0.1

pv=1-pchisq(17,12)
pv

[1] 0.1496

由于pv > 0.05，可以认为模型是可以成立的。 模型为：dec2 == 0.009