<jupyter_start><jupyter_text>Time Series and Forecasting	Assignment Course:	Time Series	and	forecasting	methods Program: MSc in Data Science (part-time)* Student: **<NAME>*** Code: **DS3616014*** Contact: **gperakis[at]aueb.gr**  Assignment Description The data you will have to analyze are in the eclass in the **excel data-assignment.xls**file. The **dependent variables** for which you will construct the models you are askedfor, are the **returns** of four investment vechicles **(Y1, Y2, Y3, Y4)** for the period **1/1991 - 12/2005**.The **independent variables** you will use in the models refer to **monthly** values/returns for the variables:* x1 = RUS-Rf* x2 = RUS (-1) -Rf (-1) lagged Russel index* x3 = MXUS-Rf* x4 = MEM- Rf* x5 = SMB* x6 = HML* x7 = MOM* x8 = SBGC-Rf* x9 = SBWG-Rf* x10 = LHY-Rf* x11 = DEFSPR* x12 = FRSI-Rf* x13 = GSCI-Rf* x14 = VIXfor the period **1/1991 - 12/2005**.Analyze dependent variables based on data for the period 1/1991 - 12/2004 **[You will not use the data for the period 1/2005 - 12/2005]**:1. Construct an appropriate **time series** model (AR, MA, ARMA).2. Develop an appropriate **regression** model  * In case of **autocorrelation** problem of **regression residuals**, correct the autocorrelation problem (using time series AR, MA, ARMA models).  * In case of **heteroscedasticity** problem of **regression residuals**, correct the    heteroskedasticity problem (using time-varying ARCH, GARCH models).3. Write the models you have found at questions (1) - (2). Assess the **goodness of fit** of these models based on the **AIC and BIC** information criteria.4. Based on the estimated models of questions (1) - (2), **construct forecasts** of the analyzed series for the period **1/2005 - 12/2005**, and **evaluate the forecasts** you have found by using two evaluation criteria:  * the mean square prediction error  * the Hit ratio (indicates the percentage of predictions that correctly evaluate the sign of the actual value of the dependent variable.)  Each student will have to analyze **two only dependent variables**  Exercise 2: Develop an appropriate regression modelThe constructed model will be of the form:$$y_{t+1} = a + b_{1}* x_{1, t} +  ... +  b_{k}* x_{k, t} + e_{t + 1}$$with $$e_{t + 1} \thicksim N(0, σ^2)$$  Load Libraries<jupyter_code># loads the library urca, which contains the command ur.df() 
# and applies Unit-Root Testing based on Dickey-Fuller test 
# install.packages(c('fGarch', 'urca','tseries','MASS'))
library(urca)
library(tseries)
library(MASS)
library(fGarch)
# Importing data into R can be carried out in various ways. Below, the command read.table is used:
data_all = read.table('/home/socital1/Desktop/george/AUEB/TimeSeries/time_series_in_r/Data-Assignment-Feb-2018.txt')
xdata_cols <- c('V5','V6', 'V7','V8','V9','V10','V11','V12','V13','V14','V15','V16','V17','V18')

xdata <- data_all[c(xdata_cols)]
head(xdata)<jupyter_output><empty_output><jupyter_text>Let create a time series object using the function “ts” from a vector - single time-series or a matrix - multivariate time-series. The data consist of the Y1 monthly returns from 1/1991 - 12/2004,<jupyter_code>y1 <- data_all$V1

j1 = ts(y1, frequency=12, start = c(1991, 1), end = c(2004, 12))
xfactors=ts(xdata, frequency=12,  start = c(1991, 1), end = c(2004, 12))

train_data = data_all[1: 168,]
test_data = data_all[169:180,]
attach(train_data)<jupyter_output><empty_output><jupyter_text>Performing Stepwise Linear Regression will all the parameters, in order to find the best model with the lower AIC<jupyter_code>fitnull <- lm(V1 ~ 1)

stepSR <- step(fitnull, 
               scope=list(lower = ~ 1,
                          upper = ~ V5 + V6 + V7 + V8 + V9 + V10 + V11 + V12 + V13 + V14 + V15 + V16 + V17 + V18),
               direction="both",
               data=train_data)

stepSR$anova
summary(stepSR)
detach(train_data)<jupyter_output><empty_output><jupyter_text>We'll now run a diagnostic check for the residuals of the fitted regression model<jupyter_code>stepSRresidualsTS <- ts(stepSR$residuals,
                        frequency=12,
                        start = c(1991,1),
                        end = c(2004, 12))
stepSRresidualsTS
par(mfrow=c(3,2))

# set up the graphics
acf(ts(stepSRresidualsTS,freq=1), 48, main="ACF of residuals")
pacf(ts(stepSRresidualsTS,freq=1), 48, main="PACF of residuals")
acf(ts(stepSRresidualsTS^2,freq=1), 48, main="ACF of squared residuals")
pacf(ts(stepSRresidualsTS^2,freq=1), 48, main="PACF of squared residuals")
qqnorm(stepSRresidualsTS,main="Normal QQplot of residuals")
qqline(stepSRresidualsTS)
plot(j1, type="l", col='red', lwd=1, main="Time Series plot of Y1", ylab="Monthly returns")<jupyter_output><empty_output><jupyter_text>Based on the residual plots presented above:* From the ACF and PACF of the residuals, it seems that the assumptions with respect to them are **not** satisfied:  * **Restricted MA(9) (1,2,8,9 lags)** and **restricted AR(8) (1,8 lags)** must be implemented in the residuals of the stepwise regression model.* From the PACF fo the squared residuals it seems that there is an **heteroskedastisity** issue.  * We'll try to fix this by implementing a ARCH(8) model in the residuals* From the Normal  QQPlot of residuals several points are off the line. So we assume non normality.From the lower right plot of the return series we observe that the volatility of the return series is not constant over time. There is evidence of volatility clustering phenomenon, since there are periods of high and low volatility.   Models Restricted AR(8) for lags: 1, 8<jupyter_code>ar8restricted <- arima(stepSRresidualsTS, 
                       order=c(8,0,0), 
                       fixed= c(NA, 0,0,0,0,0,0,NA, NA))
ar8restricted
par(mfrow=c(3,2))
ar8restricted_residualsTS <- ts(ar8restricted$residuals, frequency=12, start = c(1991, 1), end = c(2004, 12))

# set up the graphics
acf(ts(ar8restricted_residualsTS,freq=1), 48, main="ACF on residuals of restricted AR")
pacf(ts(ar8restricted_residualsTS,freq=1), 48, main="PACF on residuals of restricted AR")
acf(ts(ar8restricted_residualsTS^2,freq=1), 48, main="ACF on squared residuals of restricted AR")
pacf(ts(ar8restricted_residualsTS^2,freq=1), 48, main="PACF on squared residuals of restricted AR")
qqnorm(ar8restricted_residualsTS,main="Normal QQplot on residuals of restricted AR")
qqline(ar8restricted_residualsTS)
Box.test(ar8restricted_residualsTS, lag=12, type="Ljung")
Box.test(ar8restricted_residualsTS^2, lag=12, type="Ljung")<jupyter_output><empty_output><jupyter_text>Restricted MA(9) (1, 2, 8, 9 lags)<jupyter_code>ma9restricted <- arima(stepSRresidualsTS, 
                       order=c(0,0,9), 
                       fixed= c(NA, NA, 0, 0, 0, 0, 0, NA, NA, NA))
ma9restricted
par(mfrow=c(3,2))
ma9restricted_residualsTS <- ts(ma9restricted$residuals, frequency=12, start = c(1991, 1), end = c(2004, 12))

# set up the graphics
acf(ts(ma9restricted_residualsTS,freq=1), 48, main="ACF on residuals of restricted MA")
pacf(ts(ma9restricted_residualsTS,freq=1), 48, main="PACF on residuals of restricted MA")
acf(ts(ma9restricted_residualsTS^2,freq=1), 48, main="ACF on squared residuals of restricted MA")
pacf(ts(ma9restricted_residualsTS^2,freq=1), 48, main="PACF on squared residuals of restricted MA")
qqnorm(ma9restricted_residualsTS,main="Normal QQplot on residuals of restricted MA")
qqline(ma9restricted_residualsTS)
Box.test(ma9restricted_residualsTS, lag=12, type="Ljung")
Box.test(ma9restricted_residualsTS^2, lag=12, type="Ljung")<jupyter_output><empty_output><jupyter_text>Restricted ARMA(8, 9)<jupyter_code>arma89restricted <- arima(stepSRresidualsTS, 
                          
                       order=c(8,0,9), 
                       fixed= c(NA, 0, 0, 0, 0, 0, 0, NA, NA, NA, 0, 0,0,0, 0, NA, NA, NA))
arma89restricted
par(mfrow=c(3,2))
arma89restricted_residualsTS <- ts(arma89restricted$residuals, frequency=12, start = c(1991, 1), end = c(2004, 12))

# set up the graphics
acf(ts(arma89restricted_residualsTS,freq=1), 48, main="ACF on residuals of restricted ARMA")
pacf(ts(arma89restricted_residualsTS,freq=1), 48, main="PACF on residuals of restricted ARMA")
acf(ts(arma89restricted_residualsTS^2,freq=1), 48, main="ACF on squared residuals of restricted ARMA")
pacf(ts(arma89restricted_residualsTS^2,freq=1), 48, main="PACF on squared residuals of restricted ARMA")
qqnorm(arma89restricted_residualsTS,main="Normal QQplot on residuals of restricted MA")
qqline(arma89restricted_residualsTS)
Box.test(arma89restricted_residualsTS, lag=12, type="Ljung")
Box.test(arma89restricted_residualsTS^2, lag=12, type="Ljung")<jupyter_output><empty_output><jupyter_text>Based on the residual plots and the Box Tests of the 3 different models presented above:* From the ACF and PACF of the residuals:  * **Restricted AR**: It seems that the assumptions with respect to them are satisfied  * **Restricted MA**: It seems that the assumptions with respect to them are satisfied  * **Restricted ARMA**: It seems that the assumptions with respect to them are satisfied* From the PACF fo the **squared** residuals:  * **Restricted AR**: It seems that the assumptions regarding homoskedasticity are satisfied  * **Restricted MA**: It seems that the assumptions regarding homoskedasticity are satisfied  * **Restricted ARMA**: It seems that the assumptions regarding homoskedasticity are satisfied* From the Normal QQPlot of residuals **for all 3 models** normality is **rejected**!! Summarizing we have the following:* Restricted AR(8):   * Variance: **4.266e-05**  * AIC: **-1205.27*** Restricted MA(9):  * Variance: **4.144e-05**  * AIC: **-1205.84*** Restricted ARMA(8, 9)  * Variance: **4.065e-05**  * AIC: **-1205.05**   Conclusion:* The model with the lowest variance and lower AIC score is the best. * Here, **restricted ARMA(8, 9)** has both the smallest variance and smallest AIC. Thus, it's the model that best fits the data.  Heteroskedasticity<jupyter_code># First, we will fit an ARCH(1) model:
m1arch <- garchFit(~garch(1,0), data=stepSR$residuals, trace=F)
# trace = F reduces the summary
summary(m1arch)<jupyter_output>Title:
 GARCH Modelling 

Call:
 garchFit(formula = ~garch(1, 0), data = stepSR$residuals, trace = F) 

Mean and Variance Equation:
 data ~ garch(1, 0)
<environment: 0x6047988>
 [data = stepSR$residuals]

Conditional Distribution:
 norm 

Coefficient(s):
         mu        omega       alpha1  
-8.3325e-19   3.9679e-05   1.5778e-01  

Std. Errors:
 based on Hessian 

Error Analysis:
         Estimate  Std. Error  t value Pr(>|t|)    
mu     -8.333e-19   5.454e-04    0.000    1.000    
omega   3.968e-05   6.400e-06    6.200 5.63e-10 ***
alpha1  1.578e-01   1.240e-01    1.272    0.203    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log Likelihood:
 599.9101    normalized:  3.570894 

Description:
 Mon Mar  5 12:52:00 2018 by user:  


Standardised Residuals Tests:
                                Statistic p-Value     
 Jarque-Bera Test   R    Chi^2  4.584752  0.1010261   
 Shapiro-Wilk Test  R    W      0.9840617 0.05148097  
 Ljung-Box Test     R    Q(10)  39.710[...]<jupyter_text>We obtain the parameter estimates,the corresponding standard errors, the t-statistics and the associated p-values, in order to examine thestatistical significance of the model parameters. The ARCH(1) model is written:$$y_{t} = μ + e_{t} , $$$$e_{t} = σ_{t}*v_{t} \thicksim N(0,1), $$$$σ^2_{t} = a_{0} + a_{1}*e^2_{t-1} $$Thus, the estimated ARCH(1) model is:$$y_{t} = -1.1216e-18 + e_{t} , $$$$e_{t} = σ_{t}*v_{t} \thicksim N(0,1), $$$$σ^2_{t} = 0.000039679 + 0.15778*e^2_{t-1} $$There are different diagnostic tests outlined in the aformentioned summary, performed on the residuals $R$ and/or the squared residuals $R^2$. The **Jarque-Bera** and the **Shapiro-Wilk** tests are used in order to examine the normality assumption of the model residuals. It seems that the null hypothesis of normality (p values 0.1010261 and 0.05148097) cannot be rejected, i.e. there is evidence of **normality** of residuals. Most of the remaining tests are based on the Ljung-Box Q-statistic to examine the autocorrelation of the residuals and/or the squared residuals. The results indicate that there **exists autocorrelation** in the **residuals** series. On the contrary the tests indicate **no autocorrelation in the squared residuals**.That is, the ARCH(1) model did not capture adequately the time-varying characteristics of the return series.<jupyter_code>predict(m1arch, 6)<jupyter_output><empty_output><jupyter_text>Student t ARCH model<jupyter_code>m1arch_student=garchFit(~garch(1,0),data=stepSR$residuals, cond.dist="std",trace=F)
summary(m1arch_student)<jupyter_output>Title:
 GARCH Modelling 

Call:
 garchFit(formula = ~garch(1, 0), data = stepSR$residuals, cond.dist = "std", 
    trace = F) 

Mean and Variance Equation:
 data ~ garch(1, 0)
<environment: 0x3a8a8b0>
 [data = stepSR$residuals]

Conditional Distribution:
 std 

Coefficient(s):
         mu        omega       alpha1        shape  
-8.3325e-19   4.0900e-05   2.1345e-01   1.0000e+01  

Std. Errors:
 based on Hessian 

Error Analysis:
         Estimate  Std. Error  t value Pr(>|t|)    
mu     -8.333e-19   5.450e-04    0.000   1.0000    
omega   4.090e-05   7.825e-06    5.227 1.73e-07 ***
alpha1  2.134e-01   1.600e-01    1.334   0.1821    
shape   1.000e+01   4.108e+00    2.434   0.0149 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log Likelihood:
 598.2538    normalized:  3.561035 

Description:
 Mon Mar  5 12:52:01 2018 by user:  


Standardised Residuals Tests:
                                Statistic p-Value     
 Jarque-Bera Test   R    Chi^2  4.888336  0.086[...]<jupyter_text>The ARCH(1) model based on Student-t errors is:$$y_{t} = μ + e_{t} , $$$$e_{t} = σ_{t}*v_{t},  v_{t}\thicksim Student-t_{1.0000e+01}, $$$$σ^2_{t} = a_{0} + a_{1}*e^2_{t-1} $$Thus, the estimated ARCH(1) model is:$$y_{t} = -8.3325e-19 + e_{t} , $$$$e_{t} = σ_{t}*v_{t},  v_{t}\thicksim Student-t_{1.0000e+01}, $$$$σ^2_{t} = 0.0000409 + 0.21345*e^2_{t-1} $$There are different diagnostic tests outlined in the aformentioned summary, performed on the residuals $R$ and/or the squared residuals $R^2$. The **Jarque-Bera** and the **Shapiro-Wilk** tests are used in order to examine the normality assumption of the model residuals. It seems that the null hypothesis of normality (p values 0.08679832 and 0.03793478) cannot surely be rejected, i.e. there is *some* evidence of **normality** of residuals. Most of the remaining tests are based on the Ljung-Box Q-statistic to examine the autocorrelation of the residuals and/or the squared residuals. The results indicate that there **exists autocorrelation** in the **residuals** series. On the contrary the tests indicate some autocorrelation in the squared residuals.That is, the ARCH(1) model did not capture adequately the time-varying characteristics of the return series.  GARCH ModelsWe can estimate a GARCH(1,1) model based on normal distribution for the error process. The followingset of commands estimates the model, presents a summary of the obtained results, gives several graphsto evaluate the model fit, and computes predictions:<jupyter_code>m2garch=garchFit(~garch(1,1),data=stepSR$residuals,trace=F)
summary(m2garch)<jupyter_output>Title:
 GARCH Modelling 

Call:
 garchFit(formula = ~garch(1, 1), data = stepSR$residuals, trace = F) 

Mean and Variance Equation:
 data ~ garch(1, 1)
<environment: 0x1fdb868>
 [data = stepSR$residuals]

Conditional Distribution:
 norm 

Coefficient(s):
         mu        omega       alpha1        beta1  
-8.3325e-19   3.9679e-05   1.5778e-01   1.0000e-08  

Std. Errors:
 based on Hessian 

Error Analysis:
         Estimate  Std. Error  t value Pr(>|t|)
mu     -8.333e-19   5.584e-04    0.000    1.000
omega   3.968e-05   3.177e-05    1.249    0.212
alpha1  1.578e-01   1.704e-01    0.926    0.355
beta1   1.000e-08   7.692e-01    0.000    1.000

Log Likelihood:
 599.9101    normalized:  3.570894 

Description:
 Mon Mar  5 12:52:01 2018 by user:  


Standardised Residuals Tests:
                                Statistic p-Value     
 Jarque-Bera Test   R    Chi^2  4.584748  0.1010263   
 Shapiro-Wilk Test  R    W      0.9840617 0.05148117  
 Ljung-Box Test     R    Q(10)  39.71049  1.905[...]<jupyter_text>The GARCH(1, 1) model is written:$$y_{t} = μ + e_{t} , $$$$e_{t} = σ_{t}*v_{t},  v_{t}\thicksim N(0,1), $$$$σ^2_{t} = a_{0} + a_{1}*e^2_{t-1} +  b_{1}*σ^2_{t-1} $$Thus, the estimated GARCH(1, 1) model is:$$y_{t} = -8.3325* 10^{-19} + e_{t} , $$$$e_{t} = σ_{t}*v_{t} \thicksim N(0,1), $$$$σ^2_{t} = 3.9679 * 10^{-5} + 1.5778*10^{-1}e^2_{t-1} +  10^{-8}*σ^2_{t-1} $$There are different diagnostic tests outlined in the aformentioned summary, performed on the residuals $R$ and/or the squared residuals $R^2$. The **Jarque-Bera** and the **Shapiro-Wilk** tests are used in order to examine the normality assumption of the model residuals. It seems that the null hypothesis of normality (p values 0.1010263 and 0.05148117) cannot be rejected, i.e. there is evidence of **normality** of residuals. Most of the remaining tests are based on the Ljung-Box Q-statistic to examine the autocorrelation of the residuals and/or the squared residuals. The results indicate that there **exists autocorrelation** in the **residuals** series. On the contrary the tests indicate **no autocorrelation in the squared residuals**.That is, the GARCH(1, 1) model did not capture adequately the time-varying characteristics of the return series.  Combined ARMA + GARCH Model<jupyter_code>arma11_garch11fit <- garchFit(formula= ~arma(1,1) + garch(1,1),
                              data=stepSR$residuals,
                              trace=F)
summary(arma11_garch11fit)<jupyter_output>Title:
 GARCH Modelling 

Call:
 garchFit(formula = ~arma(1, 1) + garch(1, 1), data = stepSR$residuals, 
    trace = F) 

Mean and Variance Equation:
 data ~ arma(1, 1) + garch(1, 1)
<environment: 0x58de7a0>
 [data = stepSR$residuals]

Conditional Distribution:
 norm 

Coefficient(s):
         mu          ar1          ma1        omega       alpha1        beta1  
-8.3325e-19   8.2028e-01  -6.4688e-01   3.2707e-06   6.6956e-03   9.1483e-01  

Std. Errors:
 based on Hessian 

Error Analysis:
         Estimate  Std. Error  t value Pr(>|t|)    
mu     -8.333e-19   1.819e-04    0.000  1.00000    
ar1     8.203e-01   1.533e-01    5.350 8.78e-08 ***
ma1    -6.469e-01   1.989e-01   -3.253  0.00114 ** 
omega   3.271e-06   5.953e-06    0.549  0.58269    
alpha1  6.696e-03   3.203e-02    0.209  0.83443    
beta1   9.148e-01   1.464e-01    6.247 4.18e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log Likelihood:
 605.9441    normalized:  3.60681 

Description:
 Mon Mar[...]<jupyter_text>There are different diagnostic tests outlined in the aformentioned summary, performed on the residuals $R$ and/or the squared residuals $R^2$. The **Jarque-Bera** and the **Shapiro-Wilk** tests are used in order to examine the normality assumption of the model residuals. It seems that the null hypothesis of normality (p values 0.09175789 and 0.03454298) is rejected, i.e. there is evidence of **non normality** of residuals. Most of the remaining tests are based on the Ljung-Box Q-statistic to examine the autocorrelation of the residuals and/or the squared residuals. The results indicate that there is **no autocorrelation** in the **residuals** series. Also, the tests indicate **no autocorrelation in the squared residuals**.That is, **ARMA(1,1) + GARCH(1,1)** model **did capture** adequately the time-varying characteristics of the return series.  Forecast and Evaluation<jupyter_code>forecast <- predict(arma11_garch11fit, 12)

# plot of forecasts with 1 standard error
par(mfrow=c(1,1))
UL <- forecast$meanForecast + forecast$meanError
LL <- forecast$meanForecast - forecast$meanError

predTS <- ts(forecast$meanForecast, frequency=12, start = c(2004, 12))
UL_TS <- ts(forecast$meanForecast, frequency=12, start = c(2004, 12))
LL_TS <- ts(forecast$meanForecast, frequency=12, start = c(2004, 12))

min_x <- min(stepSRresidualsTS, LL)
max_x <- max(stepSRresidualsTS, UL)

ts.plot(stepSRresidualsTS, predTS)
lines(predTS, col="red", type="o")
lines(UL_TS, col="blue", lty="dashed")
lines(LL_TS, col="blue", lty="dashed")<jupyter_output><empty_output><jupyter_text>Mean Squared Error<jupyter_code>library(hydroGOF)
MSE <- mse(forecast$meanForecast, forecast$meanError)
MSE<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12269431060495449, 'fraction_numerical': 0.07402540516863776, 'mean_word_length': 4.062330623306233, 'pattern_counts': {'":': 0, '<': 81, '<?xml version=': 0, '>': 62, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '260431', 'license': 'MIT', 'n_tokens_mistral': 8246, 'n_tokens_neox': 7080, 'n_words': 2223, 'path': '.ipynb_checkpoints/TimeSeriesAssignment_Variable_Y1_Regression_Analysis-checkpoint.ipynb', 'repo_name': 'gperakis/time_series_in_r', 'size': '20557'}