R Notebook

rm(list=ls())

download.file("http://www.apradie.com/ec2004/salesbrands3.xlsx", "salesbrands3.xlsx", mode="wb")

trying URL 'http://www.apradie.com/ec2004/salesbrands3.xlsx'
Content type 'application/vnd.openxmlformats-officedocument.spreadsheetml.sheet' length 18194 bytes (17 KB)
==================================================
downloaded 17 KB

library(readxl)

sales_brands <- read_excel("salesbrands3.xlsx", sheet = "Resultado")

# Always familiarize with the data
head(sales_brands, 5)

# Load library xts. Install the package if you haven't already.
library(xts)

# We save the date column in a new dataset:
month <- sales_brands$date

# We delete the date column from the sales data. The date is the column 1 of the dataset: 
sales_brands <- sales_brands[,c(-1)]
# As you see, when you specify c(-1) in the column index that means that you do not want column 1! 

# Now we create an xts dataset using the previous datasets: 
sales.xts <- xts(x = sales_brands, order.by = month)

head(sales.xts, 5)

Warning: object timezone ('UTC') is different from system timezone ('')

           qbrand1tons salesbrand1 qbrand2tons salesbrand2 qbrand3tons salesbrand3
2015-01-01      106.62    60780.72      195.46    60700.51      293.10    75196.35
2015-02-01      105.82    60411.47      216.09    68289.15      494.46   105010.90
2015-03-01       88.82    53174.89      220.15    62529.54      356.88    90902.84
2015-04-01       86.69    50861.69      170.74    50575.68      291.07    72853.34
2015-05-01       87.23    51014.23      177.17    54327.48      318.09    75277.10
           qbrand4tons salesbrand4 qbrand5tons salesbrand5 qbrand6tons salesbrand6
2015-01-01      185.91    48228.51       28.94     2295.50      338.38    55561.26
2015-02-01      149.30    39922.14       24.66     1973.43      284.52    51469.74
2015-03-01      144.65    39788.22       23.75     1903.24      251.34    50395.48
2015-04-01      161.43    42692.59       21.00     1684.28      230.29    49124.61
2015-05-01      141.09    37167.52       23.44     1883.41      250.45    46906.63
           tonstotal salestotal
2015-01-01   1148.41   302762.8
2015-02-01   1274.85   327076.8
2015-03-01   1085.59   298694.2
2015-04-01    961.22   267792.2
2015-05-01    997.47   266576.4

# Do you notice any difference between data frames and xts objects?

#delete month variable
rm(month)

plot(sales.xts$qbrand2tons)

plot(sales.xts$qbrand6tons)

ln_sales <- log(sales.xts)

sales_monthly_growth_q2 = diff(ln_sales$qbrand2tons)
sales_monthly_growth_q6 = diff(ln_sales$qbrand6tons)

plot(sales_monthly_growth_q2, type = "l", col = "blue")

lines(sales_monthly_growth_q6, type = "l", col = "red")

sales_annual_growth_q2 = diff(ln_sales$qbrand2tons,lag = 12)
sales_annual_growth_q6 = diff(ln_sales$qbrand6tons,lag = 12)

plot(sales_annual_growth_q2, type = "l", col = "black")

lines(sales_annual_growth_q6, type = "l", col = "blue")

library(tseries)

adf.test(na.omit(sales_annual_growth_q2), alternative = "stationary")

    Augmented Dickey-Fuller Test

data:  na.omit(sales_annual_growth_q2)
Dickey-Fuller = -3.574, Lag order = 3, p-value = 0.04701
alternative hypothesis: stationary

#3 CHALLENGE 1 #RUN AND INTERPRET A DICKY-FULLER TEST FOR THE ANNUAL % DIFFERENCE OF BRAND 6, AND PROVIDE AN INTERPRETATION

# Load necessary library
library(tseries)

# Perform the ADF test on the annual growth rate of Brand 6
adf_test_q6 <- adf.test(na.omit(sales_annual_growth_q6), alternative = "stationary")

# Print the test result
print(adf_test_q6)

    Augmented Dickey-Fuller Test

data:  na.omit(sales_annual_growth_q6)
Dickey-Fuller = -2.3423, Lag order = 3, p-value = 0.4387
alternative hypothesis: stationary

#Con el valor de -2.3423 podemos decir que no es estacionarionaria (tiene una única raíz unitaria) y con este p-value we accept the nule hypothesis que está dentro del 95% de confianza.

library(wbstats)

# Mexico - Infant mortality
infantm<-wb_data(indicator = c("SP.DYN.IMRT.IN"), 
      country="MEX", start_date = 1980, end_date = 2020)

# Mexico - Export value
exports<-wb_data(indicator = c("TX.VAL.MRCH.XD.WD"), 
      country="MEX", start_date = 1980, end_date = 2020)

plot.ts(infantm$SP.DYN.IMRT.IN)

plot.ts(exports$TX.VAL.MRCH.XD.WD)

m1 <- lm(exports$TX.VAL.MRCH.XD.WD ~ infantm$SP.DYN.IMRT.IN)
summary(m1)

Call:
lm(formula = exports$TX.VAL.MRCH.XD.WD ~ infantm$SP.DYN.IMRT.IN)

Residuals:
    Min      1Q  Median      3Q     Max 
-21.942 -18.107  -2.227  16.959  32.253 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)            122.4937     6.8165   17.97  < 2e-16 ***
infantm$SP.DYN.IMRT.IN  -2.7632     0.2286  -12.09 9.25e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 18.19 on 39 degrees of freedom
Multiple R-squared:  0.7892,    Adjusted R-squared:  0.7838 
F-statistic:   146 on 1 and 39 DF,  p-value: 9.245e-15

#4.1 CHALLENGE 2 - QUESTION: #RESEARCH ABOUT SPURIOUS REGRESSION. With your words briefly explain in which cases you can end up with a spurious regression

When the regression model shows a high R-squared value, suggesting a good fit, which is misleading.

#When the residuals of the regression are non-stationary, indicating that the model does not properly capture the dynamics of the data. #This type of cases can ed up with a spuria regression. #When the t-statistics for the regression coefficients are significant, even though there is no true relationship between the variables.

library(tseries)
# The errors/residual of the regression are stored in the residuals
#   attribute of the m1 model:
adf.test(m1$residuals, alternative = "stationary")

    Augmented Dickey-Fuller Test

data:  m1$residuals
Dickey-Fuller = -2.5852, Lag order = 3, p-value = 0.3434
alternative hypothesis: stationary

library(quantmod)
getSymbols(Symbols=c("^MXX","^GSPC"), peridiocity = "daily",
                     from="2015-01-01", to="2017-10-02")

[1] "MXX"  "GSPC"

#Putting the adjusted prices into one dataset:
data = Ad(merge(MXX,GSPC))
data = na.omit(data)

firstmxx = as.numeric(MXX$MXX.Adjusted[1])
firstus = as.numeric(GSPC$GSPC.Adjusted[1])

invmxx = data$MXX.Adjusted / firstmxx
invus = data$GSPC.Adjusted / firstus 

inv <- merge(invus,invmxx)
plot(inv)

library(quantmod)
getSymbols(Symbols=c("^MXX","^GSPC"), peridiocity = "daily",
                     from="2017-10-03", to="2018-02-28")

[1] "MXX"  "GSPC"

#Putting the adjusted prices into one dataset:
data = Ad(merge(MXX,GSPC))
data = na.omit(data)

secondmxx = as.numeric(MXX$MXX.Adjusted[1])
secondstus = as.numeric(GSPC$GSPC.Adjusted[1])

invmxx = data$MXX.Adjusted / secondmxx
invus = data$GSPC.Adjusted / secondstus 

inv <- merge(invus,invmxx)
plot(inv)

#For both time periods, run a cointegration test and INTERPRET your results.