Homework

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

library("fBasics")

载入需要的程辑包：timeDate
载入需要的程辑包：timeSeries

1. Consider the daily stock returns of American Express (AXP), Caterpillar (CAT), and Starbucks (SBUX) from January 1999 to December 2008. The data are simple returns given in the file d- 3stocks9908.txt (date, axp, cat, sbux).

df = data.table::fread('/Users/aisling/Documents/daydayup/nus_mqf/Financial Time Series/homework/Homework/HW 1/d-3stocks9908.txt',index = 'date')
head(df)

NA

1. Compute the sample mean, standard deviation, skewness, excess kurtosis, minimum, and maximum of each simple return series. (Hint: use the R command basicStats of fBasics)

basicStats(df)

2. Transform the simple returns to log returns. Compute the sample mean, standard deviation, skewness, excess kurtosis, minimum, and maximum of each log return series.

logR = (log(df+1))
basicStats(logR)

NA

3. Test the null hypothesis that the mean of the log returns of AXP is zero. (Hint: use the R command t.test)

axp_pvalue = t.test(logR$axp)$p.value
sprintf('the p-value is %s, so we should not reject the null hypothesis that the mean of the log returns of AXP is zero.',axp_pvalue)

[1] "the p-value is 0.752367158391223, so we should not reject the null hypothesis that the mean of the log returns of AXP is zero."

4. Obtain the histogram (with nclass=40) and sample density plot of the daily log returns of AXP stock.

require(ggplot2)
ggplot(data = logR, aes(logR$axp)) +
  geom_histogram(bins = 40) +
  labs(title="sample density plot of the daily log returns of AXP stock") +
  labs(x=" log return ", y="Freq")

5. Test if log return of AXP follows normal distribution by using at least two methods. (Need to give explanation)

# 1.The Jarque-Bera Test,a type of Lagrange multiplier test, is a test for normality. It is usually used for large data sets, because other normality tests are not reliable when n is large (for example, Shapiro-Wilk isn’t reliable with n more than 2,000).A normal distribution has a skew of zero (i.e. it’s perfectly symmetrical around the mean) and a kurtosis of three; In general, a large J-B value indicates that errors are not normally distributed.

# 2. Shapiro-Wilk test is a way to tell if a random sample comes from a normal distribution. The test gives us a W value; small values indicate the sample is not normally distributed (you can reject the null hypothesis that your population is normally distributed if your values are under a certain threshold). 

axp_jb_test = normalTest(logR$axp,method="jb")
axp_jb_pvalue = axp_jb_test@test$p.value
sprintf("P Value is %s, hence %s null hypothesis of normality", axp_jb_pvalue ,ifelse(pvalue <= 0.05, "Reject","Do not reject"))

[1] "P Value is 0, hence Do not reject null hypothesis of normality"

 
axp_sw_test = normalTest(logR$axp,method="sw")
axp_sw_pvalue = axp_sw_test@test$p.value
sprintf("P Value is %s, hence %s null hypothesis of normality", axp_sw_pvalue ,ifelse(pvalue <= 0.05, "Reject","Do not reject") )

[1] "P Value is 7.78170086263107e-33, hence Do not reject null hypothesis of normality"

2. Answer the same questions as Problem 1 but using monthly returns for General Motors (GM), CRSP value-weighted index (VW), CRSP equal-weighted index (EW) and S&P composite index from January 1975 to December 2008. The returns of the indexes include dividend distributions. Data file is m-gm3dx7508.txt (date, gm, vw, ew, sp).

df_2 = data.table::fread('/Users/aisling/Documents/daydayup/nus_mqf/Financial Time Series/homework/Homework/HW 1/m-gm3dx7508.txt',index = 'date')
head(df_2)

basicStats(df_2)

NA

logR2 = (log(df_2+1))
basicStats(logR2)

# Test the null hypothesis that the mean of the log returns of stock is zero

for (col_name in colnames(logR2)[-1]){
  axp_pvalue2 = t.test(logR2[,col_name,with = F])$p.value
  res <- ifelse(axp_pvalue2 < 0.05 ,'reject','accept')
  print(sprintf('the p-value is %s, so we should %s the null hypothesis that the mean of the log returns of %s is zero.',axp_pvalue2,res,col_name))
  
  }

[1] "the p-value is 0.816604171543212, so we should accept the null hypothesis that the mean of the log returns of gm is zero."
[1] "the p-value is 7.33172261828164e-05, so we should reject the null hypothesis that the mean of the log returns of vw is zero."
[1] "the p-value is 3.42511340170795e-05, so we should reject the null hypothesis that the mean of the log returns of ew is zero."
[1] "the p-value is 0.00393916498183515, so we should reject the null hypothesis that the mean of the log returns of sp is zero."

# Obtain the histogram (with nclass=40) and sample density plot of the daily log returns of the stock.
for (col_name in colnames(logR2)[-1]){
  
  print( ggplot(data = logR2, aes(logR2[[col_name]])) +
           geom_histogram(bins = 40) +
           labs(title= cat(" density plot ")) +
           labs(x=paste0("log return of  ", col_name), y="Freq"))
  
  }

 density plot  density plot  density plot  density plot 

NA

# Test if log return of the stock follows normal distribution by using at least two methods. (Need to give explanation)

# 1.The Jarque-Bera Test,a type of Lagrange multiplier test, is a test for normality. It is usually used for large data sets, because other normality tests are not reliable when n is large (for example, Shapiro-Wilk isn’t reliable with n more than 2,000).A normal distribution has a skew of zero (i.e. it’s perfectly symmetrical around the mean) and a kurtosis of three; In general, a large J-B value indicates that errors are not normally distributed.

# 2. Shapiro-Wilk test is a way to tell if a random sample comes from a normal distribution. The test gives us a W value; small values indicate the sample is not normally distributed (you can reject the null hypothesis that your population is normally distributed if your values are under a certain threshold). 


for (col_name in colnames(logR2)[-1]){
  
  print(sprintf('============================ %s ====================================',col_name))


  jb_test = normalTest(logR2[[col_name]],method="jb")
  jb_pvalue = jb_test@test$p.value
  res_jb2 <- ifelse(axp_pvalue2 < 0.05 ,'reject','do not reject')
  print(sprintf("Jarque Bera Test: P Value is %s, hence %s null hypothesis of normality", jb_pvalue ,res_jb2))
  
  
  sw_test = normalTest(logR2[[col_name]],method="sw")
  sw_pvalue = sw_test@test$p.value
  res_sw2 <- ifelse(axp_pvalue2 < 0.05 ,'reject','do not reject')
  print(sprintf("Shapiro-Wilk test : P Value is %s, hence %s null hypothesis of normality", sw_pvalue ,res_sw2 ))
  
}

[1] "============================ gm ===================================="
[1] "Jarque Bera Test: P Value is 0, hence reject null hypothesis of normality"
[1] "Shapiro-Wilk test : P Value is 6.01277515227384e-12, hence reject null hypothesis of normality"
[1] "============================ vw ===================================="
[1] "Jarque Bera Test: P Value is 0, hence reject null hypothesis of normality"
[1] "Shapiro-Wilk test : P Value is 8.10694896883752e-11, hence reject null hypothesis of normality"
[1] "============================ ew ===================================="
[1] "Jarque Bera Test: P Value is 0, hence reject null hypothesis of normality"
[1] "Shapiro-Wilk test : P Value is 1.04817251126734e-12, hence reject null hypothesis of normality"
[1] "============================ sp ===================================="
[1] "Jarque Bera Test: P Value is 0, hence reject null hypothesis of normality"
[1] "Shapiro-Wilk test : P Value is 4.14264029673579e-09, hence reject null hypothesis of normality"

3. Consider the monthly stock returns of value-weighted index (VW) from January 1975 to December 2008 in Problem 2. Perform the tests and draw conclusions using the 5% significance level.

1. Test H0: μ = 0 versus Ha : μ ≠ 0,where μ denotes the mean return.

pvalue = t.test(logR2[['vw']] ) $p.value
sprintf("The p-value is %s, hence we %s the null hypothesis that the mean is 0",  pvalue , ifelse(pvalue <= 0.05, "Reject","Do not reject") )

[1] "The p-value is 7.33172261828164e-05, hence we Reject the null hypothesis that the mean is 0"

2. Test H0: m3 = 0 versus Ha : m3 ≠ 0, where m3 denotes the skewness.

library(moments)
# agostino.test(x)  # 峰度检验

pvalue = agostino.test(logR2[['vw']] )$p.value
sprintf("The p-value is %s, hence we %s the null hypothesis that skewness is 0", pvalue , ifelse(pvalue <= 0.05, "Reject","Do not reject") )

[1] "The p-value is 1.28785870856518e-13, hence we Reject the null hypothesis that skewness is 0"

3. Test H0: K = 3 versus Ha : K ≠ 3, where K denotes the kurtosis.

# anscombe.test(x)  # 偏度检验
pvalue = anscombe.test(logR2[['vw']] )$p.value
sprintf("The p-value is %s, hence we %s the null hypothesis that the kurtosis is 3",  pvalue ,ifelse(pvalue <= 0.05, "Reject","Do not reject") )

[1] "The p-value is 3.91860472957707e-11, hence we Reject the null hypothesis that the kurtosis is 3"

4. Consider the daily log returns of AXP stock from January 1999 to December 2008 as in Problem 1. Perform the following tests:

1. Test the null hypothesis that the skewness measure of the returns is zero;

pvalue = agostino.test(logR[['axp']] )$p.value
sprintf("The p-value is %s, hence we %s the null hypothesis that skewness is 0", pvalue , ifelse(pvalue <= 0.05, "Reject","Do not reject") )

[1] "The p-value is 1.67548197538281e-11, hence we Reject the null hypothesis that skewness is 0"

2. Test the null hypothesis that the excess kurtosis of the returns is zero.

 
pvalue = anscombe.test(logR[['axp']] )$p.value
sprintf("The p-value is %s, hence we %s the null hypothesis that the kurtosis is 3",  pvalue , ifelse(pvalue <= 0.05, "Reject","Do not reject") )

[1] "The p-value is 1.57492483442182e-76, hence we Reject the null hypothesis that the kurtosis is 3"

5. Daily foreign exchange rates (spot rates) can be obtained from the Federal Reserve Bank in St Louis (FRED). The data are the noon buying rates in New York City certified by the Federal Reserve Bank of New York. Consider the exchange rates between the U.S. dollar and the Euro from January 4, 1999 to March 8, 2013. See the file d-exuseu.txt.

1. Compute the daily log return of the exchange rate.

df_5 = data.table::fread('/Users/aisling/Documents/daydayup/nus_mqf/Financial Time Series/homework/Homework/HW 1/d-exuseu.txt')
head(df_5)

logR5 = diff(log(df_5$VALUE))
head(logR5)

[1] -0.004412021 -0.010600202  0.003089071 -0.010161114 -0.001732502  0.001213067

2. Compute the sample mean, standard deviation, skewness, excess kurtosis, minimum, and maximum of the log returns of the exchange rate.

basicStats(logR5)

3. Obtain a density plot of the daily long returns of Dollar-Euro exchange rate.

plot(density(logR5))

4. Test H0: μ = 0 versus Ha : μ ≠ 0, where μ denotes the mean of the daily log return of Dollar- Euro exchange rate.

pvalue = t.test(logR5)$p.value
sprintf("The p-value is %s, hence we %s the null hypothesis that the mean of Dollar Euro Exchange Return is 0", pvalue , ifelse(pvalue <= 0.05, "Reject","Do not reject"))

[1] "The p-value is 0.806554013958007, hence we Do not reject the null hypothesis that the mean of Dollar Euro Exchange Return is 0"

6. Consider SP500 and IBM data (in file ‘SP.csv’ and ‘IBM.csv’ respectively):

df_sp=  data.table::fread('/Users/aisling/Documents/daydayup/nus_mqf/Financial Time Series/homework/Homework/HW 1/SP.CSV',index = 'V1')
df_ibm = data.table::fread('/Users/aisling/Documents/daydayup/nus_mqf/Financial Time Series/homework/Homework/HW 1/IBM.CSV',index = 'V1')

head(df_sp)

head(df_ibm)

NA

1. Find the log return of SP and IBM by using ‘Adjusted Price’

logR_sp = diff(log(df_sp[["SP.Adjusted"]]))
logR_ibm = diff(log(df_ibm[["IBM.Adjusted"]]))

head(logR_sp)

[1] -0.052185651 -0.026072347  0.023969244  0.018511805 -0.016934050 -0.008178439

head(logR_ibm)

[1]  0.010635805 -0.009094314  0.015077732  0.011760352 -0.011861601 -0.002429755

2. Merge the above two return series into a new dataset (hint: using ‘data.frame’) and plot them (y is IBM, x is SP)

df_new = data.frame(y = logR_ibm , x = logR_sp)
qplot(df_new$y, df_new$y)

3. Find the best linear model by using AIC as criterion (need to show all your models)

m1 <- lm(y ~ x, df_new)
summary(m1)

Call:
lm(formula = y ~ x, data = df_new)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.086931 -0.006241  0.000360  0.006798  0.093286 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.0001795  0.0002345   0.766    0.444    
x           0.2314457  0.0109979  21.045   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01308 on 3111 degrees of freedom
Multiple R-squared:  0.1246,    Adjusted R-squared:  0.1243 
F-statistic: 442.9 on 1 and 3111 DF,  p-value: < 2.2e-16

sprintf('fit the model y ~ x ,AIC = %s ',AIC(m1))

[1] "fit the model y ~ x ,AIC = -18160.9543106046 "

confint(m1)

                    2.5 %       97.5 %
(Intercept) -0.0002802095 0.0006392583
x            0.2098817677 0.2530095436

# generate some other variables : 
idx = which(df_new$x <= 0)
nsp <- rep(0,length(df_new$x))
nsp[idx] = df_new$x[idx]
c1 <- rep(0,length(df_new$x))
c1[idx] = 1
 
df_new2 <- data.frame(y = df_new$y, x = df_new$x, c1, nsp)   

m2 <- lm(y ~ x + c1 , df_new2) 
summary(m2)

Call:
lm(formula = y ~ x + c1, data = df_new2)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.086854 -0.006247  0.000366  0.006818  0.093189 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 9.841e-05  4.065e-04   0.242    0.809    
x           2.341e-01  1.550e-02  15.107   <2e-16 ***
c1          1.614e-04  6.608e-04   0.244    0.807    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01308 on 3110 degrees of freedom
Multiple R-squared:  0.1246,    Adjusted R-squared:  0.1241 
F-statistic: 221.4 on 2 and 3110 DF,  p-value: < 2.2e-16

m3 <- lm(y ~ x + nsp, df_new2) 
summary(m3)

Call:
lm(formula = y ~ x + nsp, data = df_new2)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.086771 -0.006218  0.000382  0.006816  0.092673 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.566e-06  3.303e-04   0.008    0.994    
x            2.432e-01  1.900e-02  12.803   <2e-16 ***
nsp         -2.357e-02  3.099e-02  -0.760    0.447    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01308 on 3110 degrees of freedom
Multiple R-squared:  0.1248,    Adjusted R-squared:  0.1242 
F-statistic: 221.7 on 2 and 3110 DF,  p-value: < 2.2e-16

m4 <- lm(y ~ x + c1 + nsp, df_new2)
summary(m4)

Call:
lm(formula = y ~ x + c1 + nsp, data = df_new2)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.086690 -0.006189  0.000400  0.006812  0.092568 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -8.292e-05  4.707e-04  -0.176    0.860    
x            2.461e-01  2.202e-02  11.173   <2e-16 ***
c1           1.685e-04  6.609e-04   0.255    0.799    
nsp         -2.368e-02  3.100e-02  -0.764    0.445    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01308 on 3109 degrees of freedom
Multiple R-squared:  0.1248,    Adjusted R-squared:  0.124 
F-statistic: 147.8 on 3 and 3109 DF,  p-value: < 2.2e-16

library(MASS)
m5 <- stepAIC(m4) 

Start:  AIC=-26993.91
y ~ x + c1 + nsp

       Df Sum of Sq     RSS    AIC
- c1    1 0.0000111 0.53230 -26996
- nsp   1 0.0000999 0.53238 -26995
<none>              0.53228 -26994
- x     1 0.0213730 0.55366 -26873

Step:  AIC=-26995.84
y ~ x + nsp

       Df Sum of Sq     RSS    AIC
- nsp   1  0.000099 0.53239 -26997
<none>              0.53230 -26996
- x     1  0.028054 0.56035 -26838

Step:  AIC=-26997.27
y ~ x

       Df Sum of Sq     RSS    AIC
<none>              0.53239 -26997
- x     1   0.07579 0.60818 -26585

summary(m5)

Call:
lm(formula = y ~ x, data = df_new2)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.086931 -0.006241  0.000360  0.006798  0.093286 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.0001795  0.0002345   0.766    0.444    
x           0.2314457  0.0109979  21.045   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01308 on 3111 degrees of freedom
Multiple R-squared:  0.1246,    Adjusted R-squared:  0.1243 
F-statistic: 442.9 on 1 and 3111 DF,  p-value: < 2.2e-16

sprintf('fit the model y ~ x + c1 ,AIC = %s ',AIC(m2))

[1] "fit the model y ~ x + c1 ,AIC = -18159.0140536929 "

sprintf('fit the model y ~ x + nsp ,AIC = %s ',AIC(m3))

[1] "fit the model y ~ x + nsp ,AIC = -18159.5331433622 "

sprintf('fit the model y ~ x + c1 + nsp ,AIC = %s ',AIC(m4))

[1] "fit the model y ~ x + c1 + nsp ,AIC = -18157.5982207751 "

sprintf('fit the model y ~ x + c1 + nsp ,remove nonsignificant variables one at a time while minimising AIC, AIC = %s ',AIC(m5))

[1] "fit the model y ~ x + c1 + nsp ,remove nonsignificant variables one at a time while minimising AIC, AIC = -18160.9543106046 "

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.