Assignment 3

Daniel Farías | A01236327

Kathia Ruiz | A01571094

Naila Salinas | A00832702

Sofia Badillo | A02384253

Company selected: Iberdrola SA (IBDRY)

library(xts)
library(dplyr)
library(zoo)
library(tseries)
library(stats)
library(forecast)
library(astsa)
library(corrplot)
library(AER)
library(vars)
library(dynlm)
library(vars)
library(TSstudio)
library(tidyverse)
library(sarima)
library(dygraphs)

a. Background

Iberdrola SA is one of the biggest renewable energy companies in the world, with a Spain origin and 170 years of history. It is the world’s largest wind power producer and one of the world’s largest electric utilities by market capitalization. It is focus in energy transition to to reduce the negative effects of climate change and the need for a clean, reliable and smart business model. (Iberdrola SA, 2023)

Consumer Price Index is the official measure of inflation of consumer prices of the United Kingdom, based on 700 different goods and services excluding the cost of housing. (Iberdrola SA, 2023)

The variable´s performance that would be analyzed is the “Adjusted Close Price”. Stock prices usually fluctuate, they have an opening, high, low and closing price for each day. The closing price, is the price at which a trading session ends. However, there are instances when certain events get to alter the stock price, that’s when the adjusted close price is needed. It helps determine the fair value of the stock and the maintenance of an accurate historical record of the stock’s performance. (ICICI, 2023)

b. Visualization

Plot the stock price

# Import database
#file.choose()
IB <- read.csv("C:\\Users\\danyb\\OneDrive - Instituto Tecnologico y de Estudios Superiores de Monterrey\\Docs\\Documentos\\Business Intelligence\\Quinto Semestre\\Introduction to Econometrics\\ts_stock_prices\\ts_stock_prices\\IBDRY.csv")  

# setting time series format 
str(IB)

## 'data.frame':    417 obs. of  7 variables:
##  $ Date     : chr  "2015-01-05" "2015-01-12" "2015-01-19" "2015-01-26" ...
##  $ Open     : num  26.4 26.1 26.4 27.5 27.2 ...
##  $ High     : num  26.4 26.2 27.6 28.1 28.2 ...
##  $ Low      : num  25.7 25.5 26.2 27.4 26.9 ...
##  $ Close    : num  26.1 26 27.3 27.5 27 ...
##  $ Adj.Close: num  17.9 17.8 18.7 18.8 18.5 ...
##  $ Volume   : int  156400 193600 202600 300400 155400 161900 95600 161500 185400 137500 ...

It is possible to see that the variable “Date” is on character format. It is necessary for this analysis to have it in Date format.

IB$Date <- as.Date(IB$Date) #"%d/%m/%Y" this is only used when having the correct format
summary(IB$Adj.Close)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   16.96   19.96   26.19   30.45   40.60   53.90

This little summary of “Adj.Close” shows important information about the company value stocks. First, it illustrates that the minimum value for the stock in all the times is $16.96, while its maximum value is $53.90.

#Time to see for NA values.
sum(is.na(IB))  

## [1] 0

colSums(is.na(IB))  

##      Date      Open      High       Low     Close Adj.Close    Volume 
##         0         0         0         0         0         0         0

There are no missing values on the dataset.

# visualizing time series data 
# time series plot 1
plot(IB$Date,IB$Adj.Close,type="l",col="steelblue", lwd=2, xlab ="Date",ylab ="Adjusted Close Price", main = "Iberdrola SA Stock Price")

This graph shows the adjusted closing price of the stock over the years.

• In 2016, it started at around $20. It is on this year when price of the stock got its lowest value.
• It began to rise before 2018, nearly hitting $30 that year.
• There was a drop in 2020, likely due to the pandemic.
• The stock reached its highest point around 2021.
• Towards the end of 2021, there was a significant decline.
• Since then, it has been relatively stable with fewer major fluctuations.

It is pretended to have more specifically information, that is why different visualization of this data are going to be done below.

# time series plot 2
IBxts<-xts(IB$Adj.Close,order.by=IB$Date)
plot(IBxts)

In this second graph we can see the performance of the variable more clearly thanks to the grids that help to better visualize the behavior of the variable through the years, but it is on the next graph where we could see very specifically information.

# time series plot 3
dygraph(IBxts, main = "Iberdrola SA Stock Price") %>% 
  dyOptions(colors = RColorBrewer::brewer.pal(4, "Dark2")) %>%
  dyShading(from = "2017-02-6",
            to = "2022-12-22", 
            color = "#FFE6E6")

This graph illustrates a significant increase in the stock price, particularly noticeable from February 6, 2017, to December 14, 2022. The price began at $19.15 and reached a final value of $44.93 in 2022, marking an approximate increase of $25.78 USD during this period. A noticeable decline can be observed starting in March 2020, confirming our theory that this drop was due to the pandemic. However, the stock began to recover in May 2020, reaching its peak value on January 4, 2021.

# time series plot 3
dygraph(IBxts, main = "Iberdrola SA Stock Price") %>% 
  dyOptions(colors = RColorBrewer::brewer.pal(4, "Dark2")) %>%
  dyShading(from = "2022-01-1",
            to = "2022-12-22", 
            color = "#FFE6E6")

To answer the question about the variable’s performance in the last year, we’ve prepared another graph, focusing on this specific period.

• It began the year with an approximate value of $42.
• In late February, there was a decline, but it quickly recovered to its starting price.
• The peak for the year was at $45.23, in May.
• However, it hit its lowest point in October.
• By December, it closed the year almost reaching a price of $45.

Decomposing time series to identify important behaviors.

# 1) Decompose a time series
# Decomposition: observed + trend + seasonal + random
# observed: data observations 
# trend: increasing / decreasing value of data observations
# seasonality: repeating short-term cycle in time series 
# noise: random variation in time series 

IBts<-ts(IB$Adj.Close,frequency=52,start=c(2015,1))
IB_ts_decompose<-decompose(IBts)
plot(IB_ts_decompose)   

In the decomposition of the time series it is possible to get a better understanding of its behavior based on its characteristics. The observed variable part shows a similar performance to the ones previously described on the time series plots. The trend part shows a more clear positive trend than what could be seen in the plots. With this, it can be inferred that while the stock may face periodic decreases, it generally maintains a favorable positive trend.The seasonal part clearly shows a pattern of dates where the stock price is more favorable. Lastly, it is possible to see similar behaviors divided in two different periods of the random part for the variable, one before 2020, and one after it.

c. Estimation

Is time series stationary?

# Stationary Test 
adf.test(IB$Adj.Close)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  IB$Adj.Close
## Dickey-Fuller = -2.3197, Lag order = 7, p-value = 0.4424
## alternative hypothesis: stationary

### H0: Non-stationary and HA: Stationary.  
### The P-Value is higher than 0.05 so it fails to Reject the H0. Time series data is non-stationary.

The p-value obtained from the test is 0.4424, which is higher than the typical significance level of 0.05. Since the p-value is greater than 0.05, we fail to reject the null hypothesis (H0), which means that the time series data (in this case, “Adj.Close”) is considered non-stationary.

Does time series show serial autocorrelation?

# Serial Autocorrelation
# ACF plots: correlation between two periods in a time series is referred as autocorrelation function (acf)
acf(IB$Adj.Close,main="Significant Autocorrelations")  # Generally, non-stationary time series data show the presence of serial autocorrelation.  

The graph suggest that there is a significant autocorrelation in the variable of the model, giving as a result that a period have a significant impact for the next period behavior, which can become a problem because the model may not be adequately capturing the structure of the time series.

Estimating models.

In order to make different models, we are going to compare possible results using “log” and “diff” to get a stationary time series data.

plot(IB$Date,log(IB$Adj.Close), type="l",col="blue", lwd=2, xlab ="Date",ylab ="log(Stock Price)", main = "IBDRY Stock Price") #convertí variable a log y es muy similar

Transforming with log is really similar to the original one.

adf.test(log(IB$Adj.Close))

## 
##  Augmented Dickey-Fuller Test
## 
## data:  log(IB$Adj.Close)
## Dickey-Fuller = -2.1693, Lag order = 7, p-value = 0.5059
## alternative hypothesis: stationary

It is still non-stationary due the p-value of 0.5059.

plot(diff(log(IB$Adj.Close)),type="l",ylab="first order difference",main = "Diff - GE Stock Price") #diferencias, ya es estacionaria

It becomes different when using “diff” and “log”. It could be stationary.

adf.test(diff(log(IB$Adj.Close)))

## Warning in adf.test(diff(log(IB$Adj.Close))): p-value smaller than printed
## p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(log(IB$Adj.Close))
## Dickey-Fuller = -8.4297, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary

We got a p-value of 0.01 using this option, which make the time series stationary. This is one of the options that we can test.

1. ARMA With log.

Transforming the variable into log, in order to get stationary time series data. When modeling, fitted.values will contain a null that doesn’t allow to continue doing operations, so we are going to omit it. This also happens in teacher’s code. This use to be because the first value of the original data does not have anything before to be compared.

# time series modeling 
summary(IB_ARMA<-arma(log(IB$Adj.Close),order=c(1,1)))

## 
## Call:
## arma(x = log(IB$Adj.Close), order = c(1, 1))
## 
## Model:
## ARMA(1,1)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.1985876 -0.0196660  0.0005597  0.0198321  0.1085765 
## 
## Coefficient(s):
##            Estimate  Std. Error  t value Pr(>|t|)    
## ar1        0.995734    0.004385  227.062   <2e-16 ***
## ma1       -0.109437    0.049321   -2.219   0.0265 *  
## intercept  0.016535    0.014798    1.117   0.2638    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Fit:
## sigma^2 estimated as 0.001151,  Conditional Sum-of-Squares = 0.48,  AIC = -1632.55

We got a negative AIC value, which could be a good parameter to evaluate the model later. It is also interesting to see the coefficients for ar1 and ma1, the ar1 coefficient is 0.995734, which could be interpreted as a great impact in the estimating model, it is also the more significant one. The cause for this results could be the non-stationary time series data.

IB_ARMA$fitted.values <- na.omit(IB_ARMA$fitted.values)

We transformed the variable to log, but in order to get the correct values in its original format, it is necessary to apply exp.

IB_estimated_stock_price<-exp(IB_ARMA$fitted.values)

Now we can plot the result.

plot(IB_estimated_stock_price)

2. ARMA With log and diff.

It is necessary to transform the variable to log and diff because of non-stationary time series data. As we saw before, this combination can transform a non-stationary time series to a stationary one.

# time series modeling 
summary(IB_ARMA2<-arma(diff(log(IB$Adj.Close),order=c(1,1))))

## 
## Call:
## arma(x = diff(log(IB$Adj.Close), order = c(1, 1)))
## 
## Model:
## ARMA(1,1)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.1965592 -0.0195813  0.0007618  0.0205598  0.1092452 
## 
## Coefficient(s):
##             Estimate  Std. Error  t value Pr(>|t|)    
## ar1        0.7819195   0.0870985    8.977   <2e-16 ***
## ma1       -0.8722668   0.0670454  -13.010   <2e-16 ***
## intercept  0.0004729   0.0002853    1.658   0.0974 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Fit:
## sigma^2 estimated as 0.001147,  Conditional Sum-of-Squares = 0.48,  AIC = -1629.88

Although we got a higher AIC value, coefficients and significance for factors changed, but if this model is calculated with stationary data, it could be better due the risk of forecasting with non-stationary. For the next part, it will be necessary to omit null values, and transform the fitted.values of the model into its original form, reversing log and first ordering difference.

We need to convert these values into a vector “c” due the different length between fitted.values and the original ones, that we will be needing it to reverse the “diff” transformation. Having this time series as vector allows to sum them up, despite of length difference. Finally, we are going to convert them into time series again, in order to make a forecast.

IB_ARMA2$fitted.values <- na.omit(IB_ARMA2$fitted.values)
fv_exp<-c(exp(IB_ARMA2$fitted.values)) #reverting "log" operation using "exp"
IBv<-c(IB$Adj.Close) #Converting the original values into a vector.
Fixed_Values<-fv_exp+IBv #reverting "diff" operation summing the original values to the differences.

## Warning in fv_exp + IBv: longitud de objeto mayor no es múltiplo de la longitud
## de uno menor

ts_Fixed_Values <- ts(Fixed_Values, start = 1, end = length(Fixed_Values), frequency = 1) #Make it time series

Now we can plot the result.

plot(ts_Fixed_Values)

ARIMA MODEL

IB_ARIMA <- Arima(IB$Adj.Close,order=c(1,1,1))
summary(IB_ARIMA)

## Series: IB$Adj.Close 
## ARIMA(1,1,1) 
## 
## Coefficients:
##           ar1     ma1
##       -0.1545  0.0740
## s.e.   0.4457  0.4489
## 
## sigma^2 = 1.377:  log likelihood = -655.81
## AIC=1317.62   AICc=1317.68   BIC=1329.71
## 
## Training set error measures:
##                      ME     RMSE       MAE       MPE     MAPE      MASE
## Training set 0.06965504 1.169198 0.7981992 0.1753621 2.526616 0.9990932
##                      ACF1
## Training set -0.003636467

plot(IB_ARIMA$fitted)

ts_IB <- IB$Adj.Close
ets_model <- ets(ts_IB)
forecast_values <- forecast(ets_model, h = 10)
print(forecast_values)

##     Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 418       44.84566 42.89646 46.79486 41.86461 47.82671
## 419       44.84566 42.24638 47.44494 40.87040 48.82092
## 420       44.84566 41.72864 47.96268 40.07859 49.61273
## 421       44.84566 41.28505 48.40628 39.40017 50.29115
## 422       44.84566 40.89056 48.80076 38.79686 50.89446
## 423       44.84566 40.53170 49.15962 38.24802 51.44330
## 424       44.84566 40.20019 49.49113 37.74103 51.95030
## 425       44.84566 39.89054 49.80078 37.26746 52.42386
## 426       44.84566 39.59889 50.09243 36.82141 52.86991
## 427       44.84566 39.32237 50.36895 36.39852 53.29280

plot(forecast_values, main = "ETS Model")

Evaluate

# model evaluation 
# Testing serial autocorrelation in regression residuals 
# Ho: There is no serial autocorrelation 
# Ha: There is serial autocorrelation
IB_ARMA_residuals<-IB_ARMA$residuals
Box.test(IB_ARMA_residuals,lag=5,type="Ljung-Box")  

## 
##  Box-Ljung test
## 
## data:  IB_ARMA_residuals
## X-squared = 2.7633, df = 5, p-value = 0.7364

# Accept the Ho. The P-value is > 0.05 indicating that ARMA model does not show residual serial autocorrelation. 

adf.test(IB_ARMA$fitted.values) #this test was being applied to the variable with fixed values.

## 
##  Augmented Dickey-Fuller Test
## 
## data:  IB_ARMA$fitted.values
## Dickey-Fuller = -2.1065, Lag order = 7, p-value = 0.5324
## alternative hypothesis: stationary

# That's why we are always going to have a similar p-value in this test, and the previous one.
# It needs to be applied to fitted.values of the model, in order to get a real p-value.

p-value for “Ljung-Box” indicates that this model does not show residual serial autocorrelation, but we can’t ignore the fact that p-value for augmented Dickey-Fuller Test shows that H0 could be not rejected, which means that time series model is non-stationary. Forecasting with non-stationary series might bring incorrect forecasts.

IB_ARMA2_residuals<-IB_ARMA2$residuals
Box.test(IB_ARMA2_residuals,lag=5,type="Ljung-Box")  

## 
##  Box-Ljung test
## 
## data:  IB_ARMA2_residuals
## X-squared = 2.634, df = 5, p-value = 0.7562

# Accept the Ho. The P-value is > 0.05 indicating that ARMA model does not show residual serial autocorrelation. 

adf.test(IB_ARMA2$fitted.values)

## Warning in adf.test(IB_ARMA2$fitted.values): p-value smaller than printed
## p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  IB_ARMA2$fitted.values
## Dickey-Fuller = -5.3041, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary

#These values are the result of using model ARMA in the original time series, transforming it into log and diff.

For this model, we got a higher p-value for residuals, indicating that this model does not have residual serial autocorrelation, in addition, this model has a p-value for ADF test of 0.01, which reject the H0 hypothesis, meaning that we have an stationary time series.

IB_ARIMA_RESIDUALS<-IB_ARIMA$residuals
Box.test(IB_ARIMA_RESIDUALS,lag=5,type="Ljung-Box")  

## 
##  Box-Ljung test
## 
## data:  IB_ARIMA_RESIDUALS
## X-squared = 2.382, df = 5, p-value = 0.7942

# Accept the Ho. The P-value is > 0.05 indicating that ARMA model does not show residual serial autocorrelation. 

adf.test(IB_ARIMA$fitted)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  IB_ARIMA$fitted
## Dickey-Fuller = -2.2549, Lag order = 7, p-value = 0.4698
## alternative hypothesis: stationary

The first case is repeated on ARIMA model, we have an ADF p-value higher than 0.05, which means that we can not reject the H0, meaning that this model might have a non-stationary time series data.

Choosing a model.

Having the results for the diagnostic tests of the three models, we can make a decision about which one to use for generating potentially accurate forecasts. Based on these results, we have decided to use model 2, which transforms the time series into a logarithmic and differenced format. This approach helps achieve a stationary time series, reducing the margin of error when making forecasts.

Forecast

# forecasting (ARMA2)
IB_ARMA2_forecast<-forecast(ts_Fixed_Values,h=5) #h=5 for 5 periods.
IB_ARMA2_forecast

##     Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 418       45.83929 43.91120 47.76738 42.89053 48.78805
## 419       45.83929 43.26796 48.41061 41.90679 49.77178
## 420       45.83929 42.75573 48.92284 41.12340 50.55517
## 421       45.83929 42.31688 49.36169 40.45224 51.22634
## 422       45.83929 41.92664 49.75193 39.85542 51.82316

When generating a forecast with this model, we can obtain an estimate of what the stock price for the next 5 periods could be. Taking into account a 95% confidence level, these values would be as follows:

• Period 1: Price between $42.89 and $48.79
• Period 2: Price between $41.91 and $49.77
• Period 3: Price between $41.12 and $50.56
• Period 4: Price between $40.45 and $51.23
• Period 5: Price between $39.85 and $51.82

The graphs for this forecast are plotted below.

plot(IB_ARMA2_forecast)

autoplot(IB_ARMA2_forecast)

References:

Iberdrola SA (2023). Iberdrola. https://www.iberdrola.com/conocenos/nuestra-empresa