#Downloading seasonally adjusted CPI Data: 

#Step One: Go to https://fred.stlouisfed.org/series/CPIAUCSL
#Step Two: Download the CSV file. 
#Step Three Upload data into R: 
#install.packages("readr") if needed.
library(readr)
Data <- read_csv("C:/Users/Nathaniel D. Mark/Downloads/CPIAUCSL.csv")

#Generating other data to work with: 
Data$Pop_Dens[1] <- 13.8
Data$Pop_Dens[1] <- 13.75
for (i in 3:length(Data$DATE)) Data$Pop_Dens[i] <- .95*Data$Pop_Dens[i-1]+.05*Data$Pop_Dens[i-2] + .05*rnorm(1)
Data$GDP[1] <- 10000
Data$GDP[2] <- 10000
for (i in 2:length(Data$DATE)) Data$GDP[i] <- Data$GDP[i-1]+10*Data$Pop_Dens[i-1] + 100*rnorm(1)
Data$CPIAUCSL <- NULL
View(Data)

#plot(Data$DATE, Data$Pop_Dens) gives us an error because DATE has a character class.
#We want to transform it into Date class (something R can read).

Data$date <- as.Date(Data$DATE, "%d/%m/%Y")

#Then our plot looks correct: 
plot(Data$date, Data$Pop_Dens)
plot(Data$date, Data$GDP)

#b) De-trending

#First, I just save our Population Density data as a vector: 
PopD <- Data$Pop_Dens
T <- length(PopD)

#First Difference
FDPopD <- PopD[2:T] - PopD[1:(T-1)]
plot(Data$date[2:T], FDPopD )


#Annualized growth 
GrowthPopD <- 1200*(log(PopD[2:T]) - log(PopD[1:(T-1)]))
GrowthPopGDP <- 1200*(log(Data$GDP[2:T]) - log(Data$GDP[1:(T-1)]))
plot(Data$date[2:T], GrowthPopD)

#Calculating Autocorrelations
T <- length(GrowthPopD)

#VERY SLOW way: manually
autocorrelation <- numeric(20)
#autocorrelation of order 0: 
cor(GrowthPopD,GrowthPopD)
#autocorrelation of order 1: 
autocorrelation[1] <- cor(GrowthPopD[2:T],GrowthPopD[1:(T-1)])
#autocorrelation of order 4: 
autocorrelation[4] <- cor(GrowthPopD[5:T],GrowthPopD[1:(T-4)])

#Slow way: for loops

#fast way: 
#use acf function
#if needed run this code: install.packages("stats")
library(stats)
acf(GrowthPopD, lag.max = 20, plot = TRUE)

#Estimating AR or ADL models: 
#My example is an ADL(1,1) with Y variable GrowthPopD and X variable GrowthPopGDP:

#Method One: OLS

Model <- lm(GrowthPopD[2:T] ~ GrowthPopD[1:(T-1)] + GrowthPopGDP[2:T] + GrowthPopGDP[1:(T-1)])
summary(Model)

#Method Two: adl function: 
#if needed run this code: install.packages("dLagM")
library(dLagM)
Model1 <- ardlDlm(y = GrowthPopD, x = GrowthPopGDP, p = 1 , q = 1)
summary(Model1)