Homework2_Yangyichen_2024121056

第二次作业

姓名：杨逸辰 学号：2024121056 日期：20250321

2.1假设一个月度债券指数的简单回报遵循MA(1)模型：

\[R_t=a_t+0.2a_{t-1},\quad\sigma_a=0.025.\] \[a_{100}=0.01\]

(a)计算在预测起点t=100的1步和2步预测值。

(b)预测误差的标准差是多少？

(c)同时计算回报序列的滞后1和滞后2的自相关系数.

2.1

(a)计算在预测起点 t=100的1步和2步预测值

\[\hat{R}_{101}=E[R_{101}]=E[a_{101}+0.2a_{100}]=0.2\cdot0.01=0.002.\] \[\hat{R}_{102}=E[R_{102}]=E[a_{102}+0.2a_{101}]=0.2\cdot E[a_{101}]=0.\]

theta<-0.02
sigma_a<-0.025
r1_101 <- theta * 0.01
r1_102 <- 0

(b)预测误差的标准差是多少

\[\mathrm{Var}(R_{101}-\hat{R}_{101})=\mathrm{Var}(a_{101})=\sigma_a^2=0.025^2=0.000625.\] \[\begin{aligned} & \mathrm{Var}(R_{102}-\hat{R}_{102})=\mathrm{Var}(a_{102}+0.2a_{101})=\sigma_a^2+0.2^2\sigma_a^2=0.025^2(1+0.2^2)=0.000625\cdot1.04= \\ & 0.00065. \end{aligned}\]

sd1_1 <- sigma_a
sd1_2 <- sqrt(sigma_a^2 * (1 + theta^2))

(c)同时计算回报序列的滞后1和滞后2的自相关系数

\[\begin{aligned} \rho_1 & =\frac{\gamma_1}{\gamma_0}=\frac{\theta^2\sigma_a^2}{\sigma_a^2(1+\theta^2)}=\frac{0.2^2}{1+0.2^2}=\frac{0.04}{1.04}\approx0.0385. \\ \rho_2 & =0. \end{aligned}\]

rho1_1 <- theta^2 / (1 + theta^2)
rho1_2 <- 0

2.2 假设某证券的日对数回报遵循以下模型:

\[r_t=0.01+0.2r_{t-2}+a_t\]

其中at是一个均值为0、方差为0.02的高斯白噪声序列。

(a)该回报序列的均值和方差分别是多少?

(b)计算rt的1阶和2阶自相关。

(c)假设 r100 = -0.01 和r99 = 0.02。在预测起点t= 100计算1步和2步预测回 报，以及预测误差的标准差。

2.2

(a)该回报序列的均值和方差分别是多少?

\[ E[r_t] = 0.01 \] \[\sigma_r^2=(0.2)^2Var(r_{t-2})+Var(a_t)\] \[\sigma_r^2=(0.2)^2\sigma_r^2+0.02\] \[\sigma_r^2=\frac{0.02}{1-(0.2)^2}=\frac{0.02}{0.96}=0.020833\]

avg = 0.01
avg

## [1] 0.01

sigma2_r <- 0.02 / (1 - 0.2^2)
sigma2_r

## [1] 0.02083333

(b)计算rt的1阶和2阶自相关。

由于模型中rt只与rt-2相关，因此1阶自相关系数为0

\[\rho_1=0\] \[\rho_2=\frac{Cov(r_t,r_{t-2})}{\sigma_r^2}=\frac{(0.2)\sigma_r^2}{\sigma_r^2}=0.2\]

rho2_1=0
rho2_2=0.2

(c)假设 r100 = -0.01 和r99 = 0.02。在预测起点t= 100计算1步和2步预测回 报，以及预测误差的标准差。

\[r_{101}=0.01+0.2r_{99}+a_{101}=0.01+0.2\times0.02+a_{101}=0.014+a_{101}\] \[\begin{aligned} r_{102} & =0.01+0.2r_{100}+a_{102}=0.01+0.2\times(-0.01)+a_{102}=0.008+a_{102} \end{aligned}\]

r2_101 <- 0.01 + 0.2 * 0.02
r2_101

## [1] 0.014

r2_102 <- 0.01 + 0.2 * (-0.01)
r2_102

## [1] 0.008

sd2=sqrt(0.02)

2.3考虑从1948年1月到2009年3月的美国月度失业率数据(文件名m-unrate.txt)，这些数据是经过季节性调整的，并来自圣路易斯联邦储备银行。

(a)构建该序列的时间序列模型，并用该模型预测2009年4月、5月、6月和7月的失业率。

(b)此外，拟合的模型是否意味着存在商业周期?为什么?注意，可能有多个模型能很好地拟合数据，你只需要一个合适的模型。

2.3

(a)构建该序列的时间序列模型，并用该模型预测2009年4月、5月、6月和7月的失业率。

# install.packages("tidyverse")
# install.packages("forecast")
# install.packages("lubridate")
# install.packages("xts")
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── 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 errors

library(forecast)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

library(lubridate)
library(xts)

## 载入需要的程序包：zoo
## 
## 载入程序包：'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.  #
## #                                                                             #
## ###############################################################################
## 
## 载入程序包：'xts'
## 
## The following objects are masked from 'package:dplyr':
## 
##     first, last

setwd("D:/3-研究生课程/金融计量经济学/课程作业/第二次作业/")
data1 <- read.table("m-unrate.txt", header = T) #m-unrate.txt已删除2009年3月之后的数据
data1$Date <- ymd(paste(data1$Year, data1$mon, "01", sep = "-"))
head(data1)

##   Year mon dd rate       Date
## 1 1948   1  1  3.4 1948-01-01
## 2 1948   2  1  3.8 1948-02-01
## 3 1948   3  1  4.0 1948-03-01
## 4 1948   4  1  3.9 1948-04-01
## 5 1948   5  1  3.5 1948-05-01
## 6 1948   6  1  3.6 1948-06-01

unemployment_rate_ts <- ts(data1$rate, start = c(min(data1$Year), min(data1$mon)), frequency = 12)
plot(unemployment_rate_ts, main = "美国月度失业率 (1948-2009)", ylab = "失业率 (%)")

fit <- auto.arima(unemployment_rate_ts)
summary(fit)

## Series: unemployment_rate_ts 
## ARIMA(2,1,2)(0,0,2)[12] 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2     sma1     sma2
##       1.2194  -0.3574  -1.2297  0.5344  -0.2533  -0.2267
## s.e.  0.2108   0.1913   0.1930  0.1422   0.0373   0.0402
## 
## sigma^2 = 0.03723:  log likelihood = 167.66
## AIC=-321.32   AICc=-321.17   BIC=-289.13
## 
## Training set error measures:
##                       ME      RMSE       MAE        MPE     MAPE      MASE
## Training set 0.004729229 0.1920246 0.1417211 0.04719491 2.667621 0.1646808
##                     ACF1
## Training set -0.01274407

forecast_result <- forecast(fit, h = 4)
autoplot(forecast_result, main = "失业率预测 (2009年4月至7月)")

residuals <- residuals(fit)
ggAcf(residuals) + ggtitle("残差自相关图")

ggPacf(residuals) + ggtitle("残差偏自相关图")

print(forecast_result)

##          Point Forecast    Lo 80    Hi 80    Lo 95     Hi 95
## Apr 2009       8.836091 8.588821 9.083361 8.457924  9.214257
## May 2009       9.030093 8.682206 9.377980 8.498045  9.562141
## Jun 2009       9.159194 8.709234 9.609154 8.471040  9.847348
## Jul 2009       9.250688 8.689196 9.812180 8.391960 10.109416

根据运行结果选择ARIMA(2,1,2)模型，从ACF与PACF图来看，残差中存在显著的周期性成分，尤其是在Lag12和24处，可以存在大约每12个月一次的商业周期。

2.4考虑基于市场资本化分位数的NYSE/AMEX/NASDAQ Decile 1、Decile 2、Decile 9和Decile 10的月度简单回报率。数据时间跨度为1970年1月至2008年12月，数据来自CRSP。

(a)对于Decile 2和Decile10的回报序列，在5%的显著性水平下检验前12个滞后自相关的零假设。得出结论。

(b)为Decile 2的回报序列构建一个ARMA模型，进行模型检验并写出拟合的模型。

(c)使用拟合的ARMA模型生成该序列的1步至12步预测及其预测标准误差。

2.4

(a)对于Decile 2和Decile10的回报序列，在5%的显著性水平下检验前12个滞后自相关的零假设。得出结论。

library(tidyverse)
library(forecast)

setwd("D:/3-研究生课程/金融计量经济学/课程作业/第二次作业/")
data2 <- read.table("m-dec12910.txt", header = T)

dec2_returns <- data2$dec2
dec10_returns <- data2$dec10

acf_dec2 <- acf(dec2_returns, lag.max = 12, plot = FALSE)
print(acf_dec2)

## 
## Autocorrelations of series 'dec2_returns', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000  0.114 -0.057 -0.017  0.031  0.049 -0.051  0.003 -0.049 -0.039  0.027 
##     11     12 
##  0.007  0.010

acf_dec10 <- acf(dec10_returns, lag.max = 12, plot = FALSE)
print(acf_dec10)

## 
## Autocorrelations of series 'dec10_returns', 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 
##  0.061  0.130

par(mfrow = c(1, 2))
acf(dec2_returns, lag.max = 12, main = "ACF of Decile 2 Returns")
acf(dec10_returns, lag.max = 12, main = "ACF of Decile 10 Returns")

lb_dec2 <- Box.test(dec2_returns, lag = 12, type = "Ljung-Box")
print(lb_dec2)

## 
##  Box-Ljung test
## 
## data:  dec2_returns
## X-squared = 14.22, df = 12, p-value = 0.2869

lb_dec10 <- Box.test(dec10_returns, lag = 12, type = "Ljung-Box")
print(lb_dec10)

## 
##  Box-Ljung test
## 
## data:  dec10_returns
## X-squared = 41.06, df = 12, p-value = 4.789e-05

在5%的显著性水平下,dec2的p值为0.2869>0.05,不能拒绝原假设,即Decile2的回报序列前12个滞后的自相关系数没有显著证据表明存在自相关。

在5%的显著性水平下,dec10的p值为4.789e-05<0.05,拒绝原假设,即Decile10的回报序列前12个滞后的自相关系数存在自相关。

(b)为Decile 2的回报序列构建一个ARMA模型，进行模型检验并写出拟合的模型。

fit <- auto.arima(dec2_returns)

# 输出模型摘要
summary(fit)

## Series: dec2_returns 
## ARIMA(1,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1    mean
##       0.1147  0.0095
## s.e.  0.0438  0.0024
## 
## sigma^2 = 0.002353:  log likelihood = 830.29
## AIC=-1654.58   AICc=-1654.53   BIC=-1641.84
## 
## Training set error measures:
##                         ME       RMSE        MAE      MPE     MAPE      MASE
## Training set -1.671711e-05 0.04841048 0.03734009 103.8472 177.2729 0.7349234
##                    ACF1
## Training set 0.00782701

checkresiduals(fit)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,0) with non-zero mean
## Q* = 8.7243, df = 9, p-value = 0.4631
## 
## Model df: 1.   Total lags used: 10

残差随机分布，没有明显的趋势或季节性;大多数滞后项都在置信区间内，表明残差接近于白噪声;且残差近似服从正态分布，这说明ARIMA(1,0,0)模型拟合良好。

(c)使用拟合的ARMA模型生成该序列的1步至12步预测及其预测标准误差。

forecast_values <- forecast(fit, h = 12)

print(forecast_values)

##     Point Forecast       Lo 80      Hi 80       Lo 95     Hi 95
## 517    0.013824178 -0.04833693 0.07598528 -0.08124303 0.1088914
## 518    0.010007999 -0.05256038 0.07257638 -0.08568208 0.1056981
## 519    0.009570439 -0.05300328 0.07214415 -0.08612780 0.1052687
## 520    0.009520269 -0.05305352 0.07209405 -0.08617808 0.1052186
## 521    0.009514516 -0.05305927 0.07208830 -0.08618383 0.1052129
## 522    0.009513857 -0.05305993 0.07208764 -0.08618449 0.1052122
## 523    0.009513781 -0.05306000 0.07208757 -0.08618457 0.1052121
## 524    0.009513772 -0.05306001 0.07208756 -0.08618457 0.1052121
## 525    0.009513771 -0.05306001 0.07208756 -0.08618458 0.1052121
## 526    0.009513771 -0.05306001 0.07208756 -0.08618458 0.1052121
## 527    0.009513771 -0.05306001 0.07208756 -0.08618458 0.1052121
## 528    0.009513771 -0.05306001 0.07208756 -0.08618458 0.1052121

autoplot(forecast_values) +
  ggtitle("Forecast of Future Returns for Decile 2") +
  xlab("Year") + ylab("Returns")

2.5考虑1970年至2008年的IBM股票日简单回报率(文件名d-ibm3dx7008.txt)。计算IBM股票日简单回报率绝对值序列的前100个滞后自相关函数(ACF)。是否存在长程依赖性的证据?为什么?

2.5

## install.packages("tseries")
library(tseries)

setwd("D:/3-研究生课程/金融计量经济学/课程作业/第二次作业/")
data3 <- read.table("d-ibm3dx7008.txt", header = T)
abs_rtn <- abs(data3$rtn)
abs_vwretd <- abs(data3$vwretd)
abs_ewretd  <- abs(data3$ewretd)
abs_sprtrn  <- abs(data3$sprtrn)

acf(abs_rtn, lag.max = 100, main = "ACF of rtn Stock Daily Absolute Returns")

在多个滞后期下，ACF值显著不为零，这表明rtn序列存在长期依赖性.

acf(abs_vwretd, lag.max = 100, main = "ACF of vwretd Stock Daily Absolute Returns")

在多个滞后期下，ACF值显著不为零，这表明vwretd序列存在长期依赖性.

acf(abs_ewretd, lag.max = 100, main = "ACF of ewretd Stock Daily Absolute Returns")

在多个滞后期下，ACF值显著不为零，这表明ewretd序列存在长期依赖性.

acf(abs_sprtrn, lag.max = 100, main = "ACF of sprtrn Stock Daily Absolute Returns")

在多个滞后期下，ACF值显著不为零，这表明sprtrn序列存在长期依赖性.

2.6考虑美国制造业部门的电力需求数据，这些数据是经过对数处理的，表示每个月固定一天的需求，并存储在文件power6.txt中。构建该序列的时间序列模型，并用拟合的模型生成1步至24步的预测。

2.6

library(tidyverse)
library(forecast)

setwd("D:/3-研究生课程/金融计量经济学/课程作业/第二次作业/")
data4 <- read.table("power6.txt", header = FALSE) %>%
  rename(Demand = V1)
ts_data4 <- ts(data4$Demand, start = c(2025, 1), frequency = 12)
plot(ts_data4, main = "Monthly Electricity Demand of Manufacturing Sector",
     ylab = "Demand (Units)", xlab = "Year")

fit <- auto.arima(ts_data4)
summary(fit)

## Series: ts_data4 
## ARIMA(0,1,1)(1,1,0)[12] 
## 
## Coefficients:
##           ma1     sar1
##       -0.5014  -0.5313
## s.e.   0.0628   0.0541
## 
## sigma^2 = 0.0005049:  log likelihood = 596.54
## AIC=-1187.07   AICc=-1186.97   BIC=-1176.49
## 
## Training set error measures:
##                        ME       RMSE        MAE         MPE     MAPE      MASE
## Training set 0.0002055022 0.02182204 0.01714808 0.001784199 0.178798 0.4714713
##                    ACF1
## Training set 0.04213773

forecast_values <- forecast(fit, h = 24)
print(forecast_values)

##          Point Forecast     Lo 80     Hi 80     Lo 95     Hi 95
## Jan 2047       9.890679  9.861883  9.919475  9.846639  9.934719
## Feb 2047       9.870136  9.837959  9.902313  9.820925  9.919347
## Mar 2047       9.902651  9.867416  9.937886  9.848763  9.956539
## Apr 2047       9.925680  9.887632  9.963728  9.867490  9.983870
## May 2047       9.948251  9.907584  9.988918  9.886056 10.010446
## Jun 2047      10.027257  9.984129 10.070384  9.961299 10.093215
## Jul 2047      10.059916 10.014461 10.105370  9.990399 10.129433
## Aug 2047      10.088482 10.040813 10.136150 10.015579 10.161384
## Sep 2047      10.098860 10.049076 10.148644 10.022722 10.174998
## Oct 2047      10.036230  9.984417 10.088043  9.956989 10.115471
## Nov 2047       9.976174  9.922409 10.029940  9.893947 10.058402
## Dec 2047       9.954410  9.898761 10.010060  9.869302 10.039519
## Jan 2048       9.951905  9.889673 10.014136  9.856729 10.047080
## Feb 2048       9.927240  9.861533  9.992948  9.826749 10.027732
## Mar 2048       9.964582  9.895573 10.033591  9.859042 10.070122
## Apr 2048       9.979680  9.907521 10.051838  9.869322 10.090037
## May 2048      10.011085  9.935908 10.086262  9.896112 10.126059
## Jun 2048      10.092744 10.014666 10.170823  9.973333 10.212156
## Jul 2048      10.130738 10.049862 10.211614 10.007048 10.254428
## Aug 2048      10.155550 10.071969 10.239130 10.027725 10.283375
## Sep 2048      10.165674 10.079474 10.251874 10.033843 10.297505
## Oct 2048      10.089778 10.001036 10.178520  9.954059 10.225497
## Nov 2048      10.027520  9.936308 10.118733  9.888022 10.167018
## Dec 2048      10.016257  9.922638 10.109876  9.873079 10.159434

autoplot(forecast_values) +
  ggtitle("Forecast of Future Electricity Demand") +
  xlab("Year") + ylab("Demand (Units)")