Advanced Econometrics HW5

Install packages

install.packages("tidyverse", repos = "https://cloud.r-project.org")

## 
## The downloaded binary packages are in
##  /var/folders/5_/389qrkvs1sd7nkp792bslx5r0000gn/T//RtmpCOd6Rt/downloaded_packages

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.1     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors

install.packages("psych", repos = "https://cloud.r-project.org")

## 
## The downloaded binary packages are in
##  /var/folders/5_/389qrkvs1sd7nkp792bslx5r0000gn/T//RtmpCOd6Rt/downloaded_packages

library(psych)

## 
## Attaching package: 'psych'
## 
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha

install.packages("forecast", repos = "https://cloud.r-project.org")

## 
## The downloaded binary packages are in
##  /var/folders/5_/389qrkvs1sd7nkp792bslx5r0000gn/T//RtmpCOd6Rt/downloaded_packages

library(forecast)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

install.packages("tseries", repos = "https://cloud.r-project.org")

## 
## The downloaded binary packages are in
##  /var/folders/5_/389qrkvs1sd7nkp792bslx5r0000gn/T//RtmpCOd6Rt/downloaded_packages

library(tseries)
install.packages("robustbase",repos = "https://cloud.r-project.org")

## 
## The downloaded binary packages are in
##  /var/folders/5_/389qrkvs1sd7nkp792bslx5r0000gn/T//RtmpCOd6Rt/downloaded_packages

library(robustbase)
install.packages("urca",repos = "https://cloud.r-project.org")

## 
## The downloaded binary packages are in
##  /var/folders/5_/389qrkvs1sd7nkp792bslx5r0000gn/T//RtmpCOd6Rt/downloaded_packages

library(urca)
install.packages("stargazer",repos = "https://cloud.r-project.org")

## 
## The downloaded binary packages are in
##  /var/folders/5_/389qrkvs1sd7nkp792bslx5r0000gn/T//RtmpCOd6Rt/downloaded_packages

library(stargazer)

## 
## Please cite as: 
## 
##  Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2.3. https://CRAN.R-project.org/package=stargazer

install.packages("dynlm",repos = "https://cloud.r-project.org")

## 
## The downloaded binary packages are in
##  /var/folders/5_/389qrkvs1sd7nkp792bslx5r0000gn/T//RtmpCOd6Rt/downloaded_packages

library(dynlm)

## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

install.packages("forecast",repos = "https://cloud.r-project.org")

## 
## The downloaded binary packages are in
##  /var/folders/5_/389qrkvs1sd7nkp792bslx5r0000gn/T//RtmpCOd6Rt/downloaded_packages

library(forecast)
install.packages("lmtest",repos = "https://cloud.r-project.org")

## 
## The downloaded binary packages are in
##  /var/folders/5_/389qrkvs1sd7nkp792bslx5r0000gn/T//RtmpCOd6Rt/downloaded_packages

library(lmtest)

Download data

data <- read_csv("UsMacro.csv")

## Rows: 271 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (4): Year, Quarter, real_gdp_cap, fed_funds_rate
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Coverting our data into time series

Creating the GDP series over time and TS graph

GDPts <- ts(data$real_gdp_cap, start = 1954, 
             end = 2021, frequency = 4 )

Lets inspect the class of our object

class(GDPts)

## [1] "ts"

Lets create a grpah for the ts

plot(GDPts, 
     xlab = "Year",
     xlim = c(1950,2023),
     ylab = "Real GDP",
     main = "U.S. Real GDP per capita (2012 dollars)",
     sub = "N = 1954:Q3 to 2022:Q1 |Source: https://fred.stlouisfed.org/series/A939RX0Q048SBEA",
     frame = TRUE, 
     col = "pink")

plot(GDPts)
lines(GDPts)

Adding a new variable to our dataframe for Log of GDP and creating a TS

data$log_real_gdp_cap <- log(data$real_gdp_cap) 

head(data)

## # A tibble: 6 × 5
##    Year Quarter real_gdp_cap fed_funds_rate log_real_gdp_cap
##   <dbl>   <dbl>        <dbl>          <dbl>            <dbl>
## 1  1954       3        16496           1.03             9.71
## 2  1954       4        16737           0.99             9.73
## 3  1955       1        17144           1.34             9.75
## 4  1955       2        17353           1.5              9.76
## 5  1955       3        17508           1.94             9.77
## 6  1955       4        17528           2.36             9.77

Creating the grpah for the TS of log GDP

Log_GDPts <- ts(data$log_real_gdp_cap, start = 1954, 
             end = 2021, frequency = 4 )

plot(Log_GDPts, 
     xlab = "Year",
     ylab = "Log of Real GDP",
     main = " Log of U.S. Real GDP per capita (2012 dollars)",
     sub = "N = 1954:Q3 to 2022:Q1 |Source: https://fred.stlouisfed.org/series/A939RX0Q048SBEA",
     frame = TRUE, 
     col = "pink")

Creating the first difference of the log of GDP per capita using the difference of the log series over time

diff_Log_GDPts <- diff(Log_GDPts)

plot(diff_Log_GDPts,
     xlab = "Year",
     ylab = "Real GDP",
     main = "First-Difference of log of U.S. Real GDP per capita (2012 dollars)",
     sub = "N = 1954:Q3 to 2022:Q1 |Source: https://fred.stlouisfed.org/series/A939RX0Q048SBEA",
     frame = TRUE, 
     col = "pink")

length(Log_GDPts)

## [1] 269

length(diff_Log_GDPts)

## [1] 268

calculate the quarterly average of this series. Based on the frequency of the data

Lest create a new variable 100* Log GDP to have a better understanding of the effect.

data$HLOGGDP <- 100 * data$log_real_gdp_cap

head(data)

## # A tibble: 6 × 6
##    Year Quarter real_gdp_cap fed_funds_rate log_real_gdp_cap HLOGGDP
##   <dbl>   <dbl>        <dbl>          <dbl>            <dbl>   <dbl>
## 1  1954       3        16496           1.03             9.71    971.
## 2  1954       4        16737           0.99             9.73    973.
## 3  1955       1        17144           1.34             9.75    975.
## 4  1955       2        17353           1.5              9.76    976.
## 5  1955       3        17508           1.94             9.77    977.
## 6  1955       4        17528           2.36             9.77    977.

Creating the new time series for the q ave and annual growth

H_Log_GDPts <- ts(data$HLOGGDP, start = 1954, 
             end = 2021, frequency = 4 )

Creating the first difference of the 100 * log of GDP per capita using the difference of the log series over time

diff_H_Log_GDPts <- diff(H_Log_GDPts)

qave <- mean(diff_H_Log_GDPts)
qave

## [1] 0.4731123

quarterss <- 4
aavgg <- (1+qave)^quarterss-1
aavgg

## [1] 3.70916

Checking if our data for Log of GDP is stationary

Checking for time trend

Lets regress Log of GDP on a constant and on time

Lets create a new collums for year and quartert together

data$Year_q <- paste(data$Year, data$Quarter, sep = "")

For the quarter effect

lmrob1 <- lmrob(log_real_gdp_cap ~ 1 + as.numeric(Year_q), data = data)

summary(lmrob1)

## 
## Call:
## lmrob(formula = log_real_gdp_cap ~ 1 + as.numeric(Year_q), data = data)
##  \--> method = "MM"
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.20532 -0.05568  0.01293  0.03972  0.07837 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -2.850e+01  8.296e-01  -34.35   <2e-16 ***
## as.numeric(Year_q)  1.959e-03  4.190e-05   46.74   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Robust residual standard error: 0.0417 
## Multiple R-squared:  0.9781, Adjusted R-squared:  0.978 
## Convergence in 24 IRWLS iterations
## 
## Robustness weights: 
##  observation 264 is an outlier with |weight| = 0 ( < 0.00037); 
##  3 weights are ~= 1. The remaining 267 ones are summarized as
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.2734  0.7918  0.8991  0.8555  0.9674  0.9989 
## Algorithmic parameters: 
##        tuning.chi                bb        tuning.psi        refine.tol 
##         1.548e+00         5.000e-01         4.685e+00         1.000e-07 
##           rel.tol         scale.tol         solve.tol          zero.tol 
##         1.000e-07         1.000e-10         1.000e-07         1.000e-10 
##       eps.outlier             eps.x warn.limit.reject warn.limit.meanrw 
##         3.690e-04         3.678e-08         5.000e-01         5.000e-01 
##      nResample         max.it       best.r.s       k.fast.s          k.max 
##            500             50              2              1            200 
##    maxit.scale      trace.lev            mts     compute.rd fast.s.large.n 
##            200              0           1000              0           2000 
##                   psi           subsampling                   cov 
##            "bisquare"         "nonsingular"         ".vcov.avar1" 
## compute.outlier.stats 
##                  "SM" 
## seed : int(0)

stargazer(lmrob1,
          title = "Log of GDP on a constant and on time",
          dep.var.caption = "DV: Log of GDP",
          covariate.labels = c("Year+quarters"),
          notes.label = "Source:N = 1954:Q3 to 2022:Q1 |Source: https://fred.stlouisfed.org/series/A939RX0Q048SBEA",
          add.lines = TRUE,
          keep.stat = c("n","rsq"),
          digits = 2,
          type = "html",
          out = "Aeco51.html")

## 
## <table style="text-align:center"><caption><strong>Log of GDP on a constant and on time</strong></caption>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td>DV: Log of GDP</td></tr>
## <tr><td></td><td colspan="1" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td>log_real_gdp_cap</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Year+quarters</td><td>0.002<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.0000)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Constant</td><td>-28.50<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.83)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">TRUE</td><td></td></tr>
## <tr><td style="text-align:left">Observations</td><td>271</td></tr>
## <tr><td style="text-align:left">R<sup>2</sup></td><td>0.98</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Source:N = 1954:Q3 to 2022:Q1 |Source: https://fred.stlouisfed.org/series/A939RX0Q048SBEA</td><td style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>

lmrob1.1 <- lmrob(log_real_gdp_cap ~ 1 + as.numeric(Year), data = data)

stargazer(lmrob1.1,
          title = "Log of GDP on a constant and on time",
          dep.var.caption = "DV: Log of GDP",
          covariate.labels = c("Year"),
          notes.label = "Source:N = 1954:Q3 to 2022:Q1 |Source: https://fred.stlouisfed.org/series/A939RX0Q048SBEA",
          add.lines = TRUE,
          keep.stat = c("n","rsq"),
          digits = 2,
          type = "html",
          out = "Aeco5.1.html")

## 
## <table style="text-align:center"><caption><strong>Log of GDP on a constant and on time</strong></caption>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td>DV: Log of GDP</td></tr>
## <tr><td></td><td colspan="1" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td>log_real_gdp_cap</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Year</td><td>0.02<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.0004)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Constant</td><td>-28.47<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.79)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">TRUE</td><td></td></tr>
## <tr><td style="text-align:left">Observations</td><td>271</td></tr>
## <tr><td style="text-align:left">R<sup>2</sup></td><td>0.98</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Source:N = 1954:Q3 to 2022:Q1 |Source: https://fred.stlouisfed.org/series/A939RX0Q048SBEA</td><td style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>

Collect the residuals of our yearly effect and create a graph

residuals1 <- residuals(lmrob1)

min(residuals1)

## [1] -0.2053212

max(residuals1)

## [1] 0.07837234

detrended <- data.frame(YEAR = data$Year, DeTrended = residuals1)

ggplot(data = detrended, aes(x = YEAR, y = DeTrended)) +
  geom_line() + labs(title = "De-Trended Series for Log of U.S. Real GDP per capita (in 2012 dollars)", y = "Residuals Log of Real GDP", x = "Years")

Checking for unit root

dftest <- ur.df(Log_GDPts, type = "drift", lags = 1)

dftest

## 
## ############################################################### 
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test # 
## ############################################################### 
## 
## The value of the test statistic is: -1.5565 18.7873

summary(dftest)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.096388 -0.003927  0.000742  0.004688  0.073650 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  0.033994   0.019009   1.788   0.0749 .
## z.lag.1     -0.002832   0.001820  -1.556   0.1208  
## z.diff.lag   0.048600   0.061258   0.793   0.4283  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01126 on 264 degrees of freedom
## Multiple R-squared:  0.01219,    Adjusted R-squared:  0.004703 
## F-statistic: 1.628 on 2 and 264 DF,  p-value: 0.1982
## 
## 
## Value of test-statistic is: -1.5565 18.7873 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.44 -2.87 -2.57
## phi1  6.47  4.61  3.79

Since we find evidence of a unit root lets test the FD of the time series

dftest2 <- ur.df(diff_Log_GDPts, type = "drift", lags = 1)

options(scipen =999)
summary(dftest2)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.097741 -0.004020  0.000445  0.004171  0.073898 
## 
## Coefficients:
##               Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)  0.0039753  0.0007929   5.013          0.000000984 ***
## z.lag.1     -0.8621445  0.0839384 -10.271 < 0.0000000000000002 ***
## z.diff.lag  -0.0926776  0.0611012  -1.517                0.131    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01123 on 263 degrees of freedom
## Multiple R-squared:  0.4821, Adjusted R-squared:  0.4782 
## F-statistic: 122.4 on 2 and 263 DF,  p-value: < 0.00000000000000022
## 
## 
## Value of test-statistic is: -10.2712 52.7533 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.44 -2.87 -2.57
## phi1  6.47  4.61  3.79

acf(Log_GDPts)

pacf(Log_GDPts)

Using the Dickey-Fuller test for unit root on the log of the GDP series

adf.test(Log_GDPts)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  Log_GDPts
## Dickey-Fuller = -1.3226, Lag order = 6, p-value = 0.8613
## alternative hypothesis: stationary

Since the data presents a unit root we will use the DF test on the first diffrence

adf.test(diff_Log_GDPts)

## Warning in adf.test(diff_Log_GDPts): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff_Log_GDPts
## Dickey-Fuller = -6.3297, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary

Extra (Lets try the pp test) on our log model

pp.test(Log_GDPts, lshort = TRUE)

## 
##  Phillips-Perron Unit Root Test
## 
## data:  Log_GDPts
## Dickey-Fuller Z(alpha) = -5.717, Truncation lag parameter = 5, p-value
## = 0.7897
## alternative hypothesis: stationary

Lest do the pp test on the diff log

pp.test(diff_Log_GDPts, lshort = TRUE)

## Warning in pp.test(diff_Log_GDPts, lshort = TRUE): p-value smaller than printed
## p-value

## 
##  Phillips-Perron Unit Root Test
## 
## data:  diff_Log_GDPts
## Dickey-Fuller Z(alpha) = -261.25, Truncation lag parameter = 5, p-value
## = 0.01
## alternative hypothesis: stationary

Let’s do the same test but for our FED funds rate variable

Creating a time Series for FED fund rate

FEDts <- ts(data$fed_funds_rate, start = 1954, 
             end = 2021, frequency = 4 )

Lets inspect the class of our object

class(FEDts)

## [1] "ts"

Lets create a grpah for the Fed ts

plot(FEDts, 
     xlab = "Year",
     xlim = c(1950,2023),
     ylab = "FED Funds Rate",
     main = "U.S. Federal Funds Effective Rate (2012 dollars)",
     sub = "N = 1954:Q3 to 2022:Q1 |Source: https://fred.stlouisfed.org/series/A939RX0Q048SBEA",
     frame = TRUE, 
     col = "purple")

Creating the first difference of the fed fund rate

diff_FEDts <- diff(FEDts)

Checking the lengths

length(FEDts)

## [1] 269

length(diff_FEDts)

## [1] 268

Creating a graph for the first difference of the FED funds rate

plot(diff_FEDts, 
     xlab = "Year",
     xlim = c(1950,2023),
     ylab = "FED Funds Rate",
     main = "First Difference U.S. Federal Funds Effective Rate (2012 dollars)",
     sub = "N = 1954:Q3 to 2022:Q1 |Source: https://fred.stlouisfed.org/series/A939RX0Q048SBEA",
     frame = TRUE, 
     col = "purple")

Calculting the average annual growth rate and the quarterly average

qave2 <- mean(diff_FEDts)
qave2

## [1] -0.003507463

aavgg2 <- (1+qave2)^quarterss-1
aavgg2

## [1] -0.01395621

Checking if our data for FED fund rate is stationary

acf(FEDts)

pacf(FEDts)

Checking for time trend in the Fed fund rate

let’s regress the fed fund rate on a constant and on time

Let’s see the effect by quarter and years

lmrob2.1 <- lmrob(fed_funds_rate ~ 1 + as.numeric(Year_q), data = data)

stargazer(lmrob2.1,
          title = "Fed Fund Rate on a constant and on time",
          dep.var.caption = "DV: Fed Fund Rate",
          covariate.labels = c("Year + quarters"),
          notes.label = "Source:N = 1954:Q3 to 2022:Q1 |Source: https://fred.stlouisfed.org/series/A939RX0Q048SBEA",
          add.lines = TRUE,
          keep.stat = c("n","rsq"),
          digits = 2,
          type = "html",
          out = "Aeco5.2.1.html")

## 
## <table style="text-align:center"><caption><strong>Fed Fund Rate on a constant and on time</strong></caption>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td>DV: Fed Fund Rate</td></tr>
## <tr><td></td><td colspan="1" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td>fed_funds_rate</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Year + quarters</td><td>-0.01<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.001)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Constant</td><td>122.31<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(19.39)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">TRUE</td><td></td></tr>
## <tr><td style="text-align:left">Observations</td><td>271</td></tr>
## <tr><td style="text-align:left">R<sup>2</sup></td><td>0.13</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Source:N = 1954:Q3 to 2022:Q1 |Source: https://fred.stlouisfed.org/series/A939RX0Q048SBEA</td><td style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>

lmrob2 <- lmrob(fed_funds_rate ~ 1 + as.numeric(Year), data = data)

summary(lmrob2)

## 
## Call:
## lmrob(formula = fed_funds_rate ~ 1 + as.numeric(Year), data = data)
##  \--> method = "MM"
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.2594 -2.2670 -0.2223  2.1405 13.1344 
## 
## Coefficients:
##                    Estimate Std. Error t value      Pr(>|t|)    
## (Intercept)      122.311944  19.393204   6.307 0.00000000116 ***
## as.numeric(Year)  -0.059397   0.009705  -6.120 0.00000000329 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Robust residual standard error: 3.057 
## Multiple R-squared:  0.1292, Adjusted R-squared:  0.126 
## Convergence in 12 IRWLS iterations
## 
## Robustness weights: 
##  21 weights are ~= 1. The remaining 250 ones are summarized as
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.02527 0.88600 0.94570 0.89450 0.97540 0.99890 
## Algorithmic parameters: 
##        tuning.chi                bb        tuning.psi        refine.tol 
##    1.547640000000    0.500000000000    4.685061000000    0.000000100000 
##           rel.tol         scale.tol         solve.tol          zero.tol 
##    0.000000100000    0.000000000100    0.000000100000    0.000000000100 
##       eps.outlier             eps.x warn.limit.reject warn.limit.meanrw 
##    0.000369003690    0.000000003678    0.500000000000    0.500000000000 
##      nResample         max.it       best.r.s       k.fast.s          k.max 
##            500             50              2              1            200 
##    maxit.scale      trace.lev            mts     compute.rd fast.s.large.n 
##            200              0           1000              0           2000 
##                   psi           subsampling                   cov 
##            "bisquare"         "nonsingular"         ".vcov.avar1" 
## compute.outlier.stats 
##                  "SM" 
## seed : int(0)

stargazer(lmrob2,
          title = "Fed Fund Rate on a constant and on time",
          dep.var.caption = "DV: Fed Fund Rate",
          covariate.labels = c("Year"),
          notes.label = "Source:N = 1954:Q3 to 2022:Q1 |Source: https://fred.stlouisfed.org/series/A939RX0Q048SBEA",
          add.lines = TRUE,
          keep.stat = c("n","rsq"),
          digits = 2,
          type = "html",
          out = "Aeco5.2.html")

## 
## <table style="text-align:center"><caption><strong>Fed Fund Rate on a constant and on time</strong></caption>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td>DV: Fed Fund Rate</td></tr>
## <tr><td></td><td colspan="1" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td>fed_funds_rate</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Year</td><td>-0.06<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.01)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Constant</td><td>122.31<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(19.39)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">TRUE</td><td></td></tr>
## <tr><td style="text-align:left">Observations</td><td>271</td></tr>
## <tr><td style="text-align:left">R<sup>2</sup></td><td>0.13</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Source:N = 1954:Q3 to 2022:Q1 |Source: https://fred.stlouisfed.org/series/A939RX0Q048SBEA</td><td style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>

Collecting the residuals to create a graph

residuals2 <- residuals(lmrob2)

detrended2 <- data.frame(YEAR = data$Year, DeTrended = residuals2)

ggplot(data = detrended2, aes(x = YEAR, y = DeTrended)) +
  geom_line() + labs(title = "Fed Fund Rate (in 2012 dollars)", y = "Fed Fund Rate", x = "Years")

Using the Dickey-Fuller test for unit root on the FED fund rate time series including drift

dftest3 <- ur.df(FEDts, type = "drift", lags = 2)

summary(dftest3)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3604 -0.1859 -0.0322  0.2536  6.6471 
## 
## Coefficients:
##             Estimate Std. Error t value    Pr(>|t|)    
## (Intercept)  0.13066    0.08318   1.571     0.11744    
## z.lag.1     -0.02870    0.01416  -2.026     0.04373 *  
## z.diff.lag1  0.30377    0.06047   5.023 0.000000941 ***
## z.diff.lag2 -0.16119    0.06104  -2.641     0.00877 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8125 on 262 degrees of freedom
## Multiple R-squared:  0.1066, Adjusted R-squared:  0.09637 
## F-statistic: 10.42 on 3 and 262 DF,  p-value: 0.000001684
## 
## 
## Value of test-statistic is: -2.0265 2.0571 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.44 -2.87 -2.57
## phi1  6.47  4.61  3.79

testing for the first diffrerence inclunding drift

dftest4 <- ur.df(diff_FEDts, type = "drift", lags = 1)

summary(dftest4)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7082 -0.1726  0.0222  0.2642  6.4248 
## 
## Coefficients:
##              Estimate Std. Error t value             Pr(>|t|)    
## (Intercept) -0.004329   0.050111  -0.086              0.93122    
## z.lag.1     -0.887039   0.074313 -11.936 < 0.0000000000000002 ***
## z.diff.lag   0.180762   0.060624   2.982              0.00314 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8173 on 263 degrees of freedom
## Multiple R-squared:  0.3962, Adjusted R-squared:  0.3916 
## F-statistic: 86.28 on 2 and 263 DF,  p-value: < 0.00000000000000022
## 
## 
## Value of test-statistic is: -11.9365 71.2399 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.44 -2.87 -2.57
## phi1  6.47  4.61  3.79

Using the Dickey-Fuller test for unit root on the FED fund rate time series but without including drift

dftest5 <- ur.df(FEDts, type = "none", lags = 1)

summary(dftest5)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7086 -0.0934  0.0498  0.3259  6.8663 
## 
## Coefficients:
##             Estimate Std. Error t value  Pr(>|t|)    
## z.lag.1    -0.013103   0.008569  -1.529     0.127    
## z.diff.lag  0.255227   0.059391   4.297 0.0000243 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8245 on 265 degrees of freedom
## Multiple R-squared:  0.07006,    Adjusted R-squared:  0.06304 
## F-statistic: 9.983 on 2 and 265 DF,  p-value: 0.00006609
## 
## 
## Value of test-statistic is: -1.5291 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

testing for the first diffrerence without inclunding drift

dftest6 <- ur.df(diff_FEDts, type = "none", lags = 1)

summary(dftest6)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7126 -0.1769  0.0179  0.2598  6.4206 
## 
## Coefficients:
##            Estimate Std. Error t value             Pr(>|t|)    
## z.lag.1    -0.88701    0.07417 -11.959 < 0.0000000000000002 ***
## z.diff.lag  0.18074    0.06051   2.987              0.00308 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8157 on 264 degrees of freedom
## Multiple R-squared:  0.3962, Adjusted R-squared:  0.3916 
## F-statistic: 86.61 on 2 and 264 DF,  p-value: < 0.00000000000000022
## 
## 
## Value of test-statistic is: -11.9587 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

Extra testing for pp test

pp.test(FEDts, lshort = TRUE)

## 
##  Phillips-Perron Unit Root Test
## 
## data:  FEDts
## Dickey-Fuller Z(alpha) = -13.458, Truncation lag parameter = 5, p-value
## = 0.3553
## alternative hypothesis: stationary

pp.test(diff_FEDts, lshort = TRUE)

## Warning in pp.test(diff_FEDts, lshort = TRUE): p-value smaller than printed
## p-value

## 
##  Phillips-Perron Unit Root Test
## 
## data:  diff_FEDts
## Dickey-Fuller Z(alpha) = -194.79, Truncation lag parameter = 5, p-value
## = 0.01
## alternative hypothesis: stationary

Calculating the dynamic effects of the Fed Fund Rate on the Log of of GDP

DEGDP <- dynlm(diff_Log_GDPts ~ diff_FEDts)

summary(DEGDP)

## 
## Time series regression with "ts" data:
## Start = 1954(2), End = 2021(1)
## 
## Call:
## dynlm(formula = diff_Log_GDPts ~ diff_FEDts)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.094196 -0.003735  0.000425  0.004475  0.067454 
## 
## Coefficients:
##              Estimate Std. Error t value         Pr(>|t|)    
## (Intercept) 0.0047447  0.0006602   7.186 0.00000000000672 ***
## diff_FEDts  0.0038755  0.0007766   4.991 0.00000108948467 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01081 on 266 degrees of freedom
## Multiple R-squared:  0.08561,    Adjusted R-squared:  0.08218 
## F-statistic: 24.91 on 1 and 266 DF,  p-value: 0.000001089

BIC <- function(DEGDP) {
  
  ssr <- sum(DEGDP$residuals^2)
  t <- length(DEGDP$residuals)
  npar <- length(DEGDP$coef)
  
  return(
    round(c("p" = npar - 1,
          "BIC" = log(ssr/t) + npar * log(t)/t,
          "R2" = summary(DEGDP)$r.squared), 4)
  )
}

BIC(dynlm(ts(diff_Log_GDPts) ~ 1))

##       p     BIC      R2 
##  0.0000 -8.9519  0.0000

order <- 1:12

BICs <- sapply(order, function(x) 
        "AR" = BIC(dynlm(ts(diff_Log_GDPts) ~ L(ts(diff_Log_GDPts), 1:x))))

BICs

##        [,1]    [,2]    [,3]    [,4]    [,5]    [,6]    [,7]    [,8]    [,9]
## p    1.0000  2.0000  3.0000  4.0000  5.0000  6.0000  7.0000  8.0000  9.0000
## BIC -8.9332 -8.9276 -8.9036 -8.8806 -8.8593 -8.8391 -8.8165 -8.7982 -8.7764
## R2   0.0031  0.0112  0.0106  0.0122  0.0158  0.0163  0.0192  0.0234  0.0262
##       [,10]   [,11]   [,12]
## p   10.0000 11.0000 12.0000
## BIC -8.7528 -8.7317 -8.7125
## R2   0.0283  0.0297  0.0367

BICs2 <- sapply(order, function(x) 
         BIC(dynlm(diff_Log_GDPts ~ L(diff_Log_GDPts, 1:x) + L(diff_FEDts, 1:x), 
                   start = c(1954, 3), end = c(2022, 1))))

BICs2

##        [,1]    [,2]    [,3]    [,4]    [,5]    [,6]    [,7]    [,8]    [,9]
## p    2.0000  4.0000  6.0000  8.0000 10.0000 12.0000 14.0000 16.0000 18.0000
## BIC -8.9144 -8.9233 -8.8814 -8.8542 -8.8204 -8.7942 -8.7505 -8.7111 -8.6769
## R2   0.0053  0.0478  0.0503  0.0679  0.0796  0.0943  0.0975  0.1021  0.1131
##       [,10]   [,11]   [,12]
## p   20.0000 22.0000 24.0000
## BIC -8.6398 -8.5977 -8.5507
## R2   0.1227  0.1252  0.1267

GDPDLM <- dynlm(diff_Log_GDPts ~ diff_FEDts + L(diff_FEDts, 2))

coeftest(GDPDLM)

## 
## t test of coefficients:
## 
##                     Estimate  Std. Error t value         Pr(>|t|)    
## (Intercept)       0.00463388  0.00065743  7.0485 0.00000000001588 ***
## diff_FEDts        0.00368513  0.00077511  4.7543 0.00000328285008 ***
## L(diff_FEDts, 2) -0.00147840  0.00077486 -1.9080          0.05749 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(GDPDLM)

## 
## Time series regression with "ts" data:
## Start = 1954(4), End = 2021(1)
## 
## Call:
## dynlm(formula = diff_Log_GDPts ~ diff_FEDts + L(diff_FEDts, 2))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.095127 -0.003742  0.000405  0.004432  0.067009 
## 
## Coefficients:
##                    Estimate Std. Error t value        Pr(>|t|)    
## (Intercept)       0.0046339  0.0006574   7.049 0.0000000000159 ***
## diff_FEDts        0.0036851  0.0007751   4.754 0.0000032828501 ***
## L(diff_FEDts, 2) -0.0014784  0.0007749  -1.908          0.0575 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01072 on 263 degrees of freedom
## Multiple R-squared:  0.09786,    Adjusted R-squared:  0.091 
## F-statistic: 14.26 on 2 and 263 DF,  p-value: 0.000001313

Let’s caluclate a point forecast for the log of the GDP series for 2022-Q2

myforecast <- forecast(Log_GDPts, level=c(95), h= 4*4)

myforecast

##         Point Forecast    Lo 95    Hi 95
## 2021 Q2       10.98350 10.96022 11.00677
## 2021 Q3       10.98819 10.95560 11.02077
## 2021 Q4       10.99287 10.95310 11.03265
## 2022 Q1       10.99756 10.95170 11.04343
## 2022 Q2       11.00225 10.95102 11.05349
## 2022 Q3       11.00694 10.95084 11.06304
## 2022 Q4       11.01163 10.95105 11.07221
## 2023 Q1       11.01632 10.95157 11.08108
## 2023 Q2       11.02101 10.95233 11.08969
## 2023 Q3       11.02570 10.95330 11.09810
## 2023 Q4       11.03039 10.95445 11.10633
## 2024 Q1       11.03508 10.95576 11.11440
## 2024 Q2       11.03977 10.95720 11.12234
## 2024 Q3       11.04446 10.95876 11.13016
## 2024 Q4       11.04915 10.96043 11.13787
## 2025 Q1       11.05384 10.96219 11.14549

plot(myforecast)

Using the arima function to create the forecast

mymodel <- auto.arima(Log_GDPts)

auto.arima(Log_GDPts, ic = "aic", trace = TRUE)

## 
##  Fitting models using approximations to speed things up...
## 
##  ARIMA(2,1,2)(1,0,1)[4] with drift         : -1624.746
##  ARIMA(0,1,0)           with drift         : -1631.35
##  ARIMA(1,1,0)(1,0,0)[4] with drift         : -1627.591
##  ARIMA(0,1,1)(0,0,1)[4] with drift         : -1628.379
##  ARIMA(0,1,0)                              : -1589.787
##  ARIMA(0,1,0)(1,0,0)[4] with drift         : -1629.995
##  ARIMA(0,1,0)(0,0,1)[4] with drift         : -1629.72
##  ARIMA(0,1,0)(1,0,1)[4] with drift         : -1632.221
##  ARIMA(0,1,0)(2,0,1)[4] with drift         : -1626.228
##  ARIMA(0,1,0)(1,0,2)[4] with drift         : -1631.875
##  ARIMA(0,1,0)(0,0,2)[4] with drift         : -1629.225
##  ARIMA(0,1,0)(2,0,0)[4] with drift         : -1627.496
##  ARIMA(0,1,0)(2,0,2)[4] with drift         : -1625.447
##  ARIMA(1,1,0)(1,0,1)[4] with drift         : -1629.138
##  ARIMA(0,1,1)(1,0,1)[4] with drift         : -1630.539
##  ARIMA(1,1,1)(1,0,1)[4] with drift         : -1628.16
##  ARIMA(0,1,0)(1,0,1)[4]                    : -1605.03
## 
##  Now re-fitting the best model(s) without approximations...
## 
##  ARIMA(0,1,0)(1,0,1)[4] with drift         : Inf
##  ARIMA(0,1,0)(1,0,2)[4] with drift         : Inf
##  ARIMA(0,1,0)           with drift         : -1640.161
## 
##  Best model: ARIMA(0,1,0)           with drift

## Series: Log_GDPts 
## ARIMA(0,1,0) with drift 
## 
## Coefficients:
##        drift
##       0.0047
## s.e.  0.0007
## 
## sigma^2 = 0.0001276:  log likelihood = 822.08
## AIC=-1640.16   AICc=-1640.12   BIC=-1632.98

mymodel

## Series: Log_GDPts 
## ARIMA(0,1,0) with drift 
## 
## Coefficients:
##        drift
##       0.0047
## s.e.  0.0007
## 
## sigma^2 = 0.0001276:  log likelihood = 822.08
## AIC=-1640.16   AICc=-1640.12   BIC=-1632.98

plot.ts(mymodel$residuals)

?forecast

## Help on topic 'forecast' was found in the following packages:
## 
##   Package               Library
##   generics              /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library
##   forecast              /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library
## 
## 
## Using the first match ...

myforecast2 <- forecast(mymodel, level = c(95), h=4*4)

plot(myforecast2,
     xlab = "Year",
     ylab = "Log of Real GDP",
     main = "Forecast for the Log of real GDP per capita",
     frame = TRUE,
     sub = "Source:N = 1954:Q3 to 2022:Q1 |Source: https://fred.stlouisfed.org/series/A939RX0Q048SBEA",
     col = "pink")

myforecast2

##         Point Forecast    Lo 95    Hi 95
## 2021 Q2       10.98355 10.96140 11.00569
## 2021 Q3       10.98828 10.95696 11.01959
## 2021 Q4       10.99301 10.95465 11.03136
## 2022 Q1       10.99774 10.95345 11.04203
## 2022 Q2       11.00247 10.95296 11.05198
## 2022 Q3       11.00720 10.95296 11.06144
## 2022 Q4       11.01193 10.95335 11.07052
## 2023 Q1       11.01666 10.95403 11.07929
## 2023 Q2       11.02139 10.95496 11.08782
## 2023 Q3       11.02613 10.95610 11.09615
## 2023 Q4       11.03086 10.95742 11.10430
## 2024 Q1       11.03559 10.95888 11.11229
## 2024 Q2       11.04032 10.96048 11.12016
## 2024 Q3       11.04505 10.96220 11.12790
## 2024 Q4       11.04978 10.96402 11.13554
## 2025 Q1       11.05451 10.96594 11.14309

dtmyforecast <- data.frame(myforecast2)

dtmyforecast$RealPoint.Forecast <- exp(dtmyforecast$Point.Forecast)

dtmyforecast

##         Point.Forecast    Lo.95    Hi.95 RealPoint.Forecast
## 2021 Q2       10.98355 10.96140 11.00569           58896.99
## 2021 Q3       10.98828 10.95696 11.01959           59176.30
## 2021 Q4       10.99301 10.95465 11.03136           59456.93
## 2022 Q1       10.99774 10.95345 11.04203           59738.90
## 2022 Q2       11.00247 10.95296 11.05198           60022.20
## 2022 Q3       11.00720 10.95296 11.06144           60306.85
## 2022 Q4       11.01193 10.95335 11.07052           60592.84
## 2023 Q1       11.01666 10.95403 11.07929           60880.19
## 2023 Q2       11.02139 10.95496 11.08782           61168.91
## 2023 Q3       11.02613 10.95610 11.09615           61458.99
## 2023 Q4       11.03086 10.95742 11.10430           61750.45
## 2024 Q1       11.03559 10.95888 11.11229           62043.29
## 2024 Q2       11.04032 10.96048 11.12016           62337.52
## 2024 Q3       11.04505 10.96220 11.12790           62633.14
## 2024 Q4       11.04978 10.96402 11.13554           62930.17
## 2025 Q1       11.05451 10.96594 11.14309           63228.61

Box.test(mymodel$resid, lag = 2, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  mymodel$resid
## X-squared = 3.6664, df = 2, p-value = 0.1599

Testing for co-integration

ur.ers(window(GDPts, c(1954), c(2022)) - window(FEDts, c(1954), c(2022)),
       model = "constant",
       lag.max = 3)

## Warning in window.default(x, ...): 'end' value not changed

## Warning in window.default(x, ...): 'end' value not changed

## 
## ############################################################### 
## # Elliot, Rothenberg and Stock Unit Root / Cointegration Test # 
## ############################################################### 
## 
## The value of the test statistic is: 3.7621

coi <- ur.df(window(GDPts, c(1954,3), c(2022,1)) - window(FEDts, c(1954,3), c(2022,1)), 
      lags = 3,
      selectlags = "AIC",
      type = "trend")

## Warning in window.default(x, ...): 'end' value not changed

## Warning in window.default(x, ...): 'end' value not changed

coi

## 
## ############################################################### 
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test # 
## ############################################################### 
## 
## The value of the test statistic is: -3.2298 13.9531 5.5151

plot(coi)

Co integration test with only drift

ur.df(window(GDPts, c(1954,3), c(2022,1)) - window(FEDts, c(1954,3), c(2022,1)), 
      lags = 3,
      selectlags = "AIC",
      type = "drift")

## Warning in window.default(x, ...): 'end' value not changed

## Warning in window.default(x, ...): 'end' value not changed

## 
## ############################################################### 
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test # 
## ############################################################### 
## 
## The value of the test statistic is: 0.4307 14.9454

Co integration test with none

ur.df(window(GDPts, c(1954,3), c(2022,1)) - window(FEDts, c(1954,3), c(2022,1)), 
      lags = 3,
      selectlags = "AIC",
      type = "none")

## Warning in window.default(x, ...): 'end' value not changed

## Warning in window.default(x, ...): 'end' value not changed

## 
## ############################################################### 
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test # 
## ############################################################### 
## 
## The value of the test statistic is: 5.2536

Extra extra testing

mymodel2 <- auto.arima(FEDts, trace = TRUE)

## 
##  Fitting models using approximations to speed things up...
## 
##  ARIMA(2,1,2)(1,0,1)[4] with drift         : Inf
##  ARIMA(0,1,0)           with drift         : 677.9387
##  ARIMA(1,1,0)(1,0,0)[4] with drift         : 669.0115
##  ARIMA(0,1,1)(0,0,1)[4] with drift         : 655.1837
##  ARIMA(0,1,0)                              : 675.9131
##  ARIMA(0,1,1)           with drift         : 653.1344
##  ARIMA(0,1,1)(1,0,0)[4] with drift         : 658.9094
##  ARIMA(0,1,1)(1,0,1)[4] with drift         : 659.086
##  ARIMA(1,1,1)           with drift         : 649.0388
##  ARIMA(1,1,1)(1,0,0)[4] with drift         : 653.3578
##  ARIMA(1,1,1)(0,0,1)[4] with drift         : 649.948
##  ARIMA(1,1,1)(1,0,1)[4] with drift         : 655.4051
##  ARIMA(1,1,0)           with drift         : 663.8734
##  ARIMA(2,1,1)           with drift         : 651.6098
##  ARIMA(1,1,2)           with drift         : 650.8487
##  ARIMA(0,1,2)           with drift         : 648.3055
##  ARIMA(0,1,2)(1,0,0)[4] with drift         : 653.7581
##  ARIMA(0,1,2)(0,0,1)[4] with drift         : 650.0828
##  ARIMA(0,1,2)(1,0,1)[4] with drift         : Inf
##  ARIMA(0,1,3)           with drift         : 649.5762
##  ARIMA(1,1,3)           with drift         : 651.4875
##  ARIMA(0,1,2)                              : 646.2487
##  ARIMA(0,1,2)(1,0,0)[4]                    : 651.6926
##  ARIMA(0,1,2)(0,0,1)[4]                    : 648.0099
##  ARIMA(0,1,2)(1,0,1)[4]                    : Inf
##  ARIMA(0,1,1)                              : 651.0919
##  ARIMA(1,1,2)                              : 648.7733
##  ARIMA(0,1,3)                              : 647.5031
##  ARIMA(1,1,1)                              : 646.9785
##  ARIMA(1,1,3)                              : 649.3971
## 
##  Now re-fitting the best model(s) without approximations...
## 
##  ARIMA(0,1,2)                              : 646.1291
## 
##  Best model: ARIMA(0,1,2)

mymodel2

## Series: FEDts 
## ARIMA(0,1,2) 
## 
## Coefficients:
##          ma1      ma2
##       0.3471  -0.1509
## s.e.  0.0587   0.0559
## 
## sigma^2 = 0.6421:  log likelihood = -320.02
## AIC=646.04   AICc=646.13   BIC=656.81

Using another package

adf.test(FEDts, alternative = "stationary",  k=5 )

## 
##  Augmented Dickey-Fuller Test
## 
## data:  FEDts
## Dickey-Fuller = -3.2444, Lag order = 5, p-value = 0.08089
## alternative hypothesis: stationary

Testing

dpLog_GDPTS <- decompose(Log_GDPts)

dpLog_GDPTS$seasonal

##                Qtr1           Qtr2           Qtr3           Qtr4
## 1954  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1955  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1956  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1957  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1958  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1959  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1960  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1961  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1962  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1963  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1964  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1965  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1966  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1967  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1968  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1969  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1970  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1971  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1972  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1973  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1974  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1975  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1976  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1977  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1978  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1979  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1980  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1981  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1982  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1983  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1984  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1985  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1986  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1987  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1988  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1989  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1990  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1991  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1992  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1993  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1994  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1995  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1996  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1997  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1998  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 1999  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 2000  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 2001  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 2002  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 2003  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 2004  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 2005  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 2006  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 2007  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 2008  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 2009  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 2010  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 2011  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 2012  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 2013  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 2014  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 2015  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 2016  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 2017  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 2018  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 2019  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 2020  0.00077239844 -0.00033228035  0.00005230922 -0.00049242731
## 2021  0.00077239844

plot(dpLog_GDPTS)

Creating a De-trended series

LogGDPdetrended <- Log_GDPts - dpLog_GDPTS$trend - dpLog_GDPTS$seasonal

plot(LogGDPdetrended)