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

#Now run a regression using these series. Report the result of the regression. Did you find significant relationship? is your result what you expected? #I don´t see any significant relation because while the infant mortality graph begins up and with a significant curve, the export value goes from a low number to a higher an is increasing while the othe is decreasing. And this is the result that i expect.

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)

#2nd period
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.

m2 <- lm( invmxx ~ invus )
summary(m2)

Call:
lm(formula = invmxx ~ invus)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.037454 -0.014008 -0.002567  0.017749  0.040777 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.79716    0.05928  13.447  < 2e-16 ***
invus        0.16206    0.05640   2.874  0.00501 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.02006 on 95 degrees of freedom
Multiple R-squared:  0.07997,   Adjusted R-squared:  0.07028 
F-statistic: 8.257 on 1 and 95 DF,  p-value: 0.005008

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

    Augmented Dickey-Fuller Test

data:  m2$residuals
Dickey-Fuller = -1.8913, Lag order = 4, p-value = 0.6214
alternative hypothesis: stationary

#Como el valor de p es mayor a 0.05 podemos decir que la hipotesis nula es aceptada con 95% de confianza que la cointegración no es estacionaria ya que no hay suficiente evidencia para demostrar lo contrario. Por el caso del análisis dickey fuller el valor del estadístico es menos negativo que el valor crítico correspondiente al nivel de significancia ( -1.8913 para el 5% de significancia), no se rechaza la hipótesis nula por que no hay evidencia suficiente. En conclusión podemos decir ya que no es estacionario que no tienen una cointegración.