R Log and Other Materials on Seasonal Adjustment for STAT 730
=============================================================
Eric Slud                                                5/5/26

> library(astsa)
  library(seasonal)


> plot(unemp)         ### length 323 monthly data, 26 years minus 1 month
> plot(log(unemp))
> plot(unemp)
> m1= seas(unemp)
> summary(m1)

Call:
seas(x = unemp)

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
AR-Nonseasonal-01  0.94360    0.03441   27.43   <2e-16 ***
MA-Nonseasonal-01  0.82540    0.05654   14.60   <2e-16 ***
MA-Seasonal-12     0.85071    0.03362   25.30   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

SEATS adj.  ARIMA: (1 1 1)(0 1 1)  Obs.: 323  Transform: none
AICc:  4324, BIC:  4339  QS (no seasonality in final):    0  
Box-Ljung (no autocorr.): 22.04   Shapiro (normality): 0.9946  

### NB models fitted:
SEATS adj.  ARIMA: (1 1 1)(0 1 1) says ARIMA(1,1,1) for nonseasonal part, (0,1,1) seasonal

> dimnames(m1$data)[[2]]
[1] "final"       "seasonal"    "seasonaladj" "trend"       "irregular"  
[6] "adjustfac"  
> 
> plot(unemp, main="Monthly Unemployment 1990-2015, with Trend")
  lines(m1$data[,4], col="blue")
  legend(locator(), legend=c("original","fitted trend"), lty=1, col=c("black","blue"))
## saved as UnempPic1.pdf

#--------------------------- 
> tmp1 = unemp - m1$data[,4] - m1$data[,5]   ### subtract "trend" and "irregular" components
  plot(tmp1, m1$data[,2])           ## exactly gives "seasonal" component

#---------------------------

> unemp.arr = array(c(unemp,NA), c(12,27))
  plot(c(unemp,NA) - rep( apply(unemp.arr,2,mean, na.rm=T), rep(12,27)), type="b")
          ### Show oscillation by removing yearly means, but amplitudes do change ...

  lgunemp.arr = array(c(log(unemp),NA), c(12,27))
  plot(1990+(0:323)/12, c(log(unemp),NA) - rep( apply(lgunemp.arr,2,mean, na.rm=T), rep(12,27)), 
       xlab = "Date", ylab="Within-Year Oscillation", type="b")
          ### same thing with log(unemp): a little more stable
          ### saved as   UnempPic2.pdf
  
### Now look at average across years to find within-year seasonal pattern

 osc =  c(log(unemp),0) - rep( apply(lgunemp.arr,2,mean, na.rm=T), rep(12,27) )
 osc[324] = NA

  plot(apply(array(osc,c(12,27)), 1, mean, na.rm=T), type="b",
       xlab="Month", ylab = "Within-year Oscillation")
          ### saved as   UnempPic3.pdf

##--------------------------------------------------------------------
### Now look at Trend as compared with moving-window-averaged filter

tmp2 = numeric(323-12)
for(i in 7:317) tmp2[i-6] = sum(c(0.5, rep(1,11), 0.5)/12 * unemp[(i-6):(i+6)])

> plot(tmp2,m1$data[7:317,"trend"])       ### close to 45-degree straight line !!
  abline(0,1, col="blue")

  plot(1990+(0:322)/12,m1$data[,"trend"], type="b", xlab="Date", 
     ylab="Unemp", main="Unemployment Trend, 1990-2017, done two ways")
  lines((1990+(0:322)/12)[7:317],tmp2, col="blue")
  legend(locator(), legend=c("Trend from `seasonal' package",
     "From simple moving average"), lty=c(NA,1), pch=c(1,NA), col=c("black","blue"))
                     ### saved as   UnempPic4.pdf

##-------------------------------------------------------------------
## WE HAVE SEEN THAT "IRREGULAR" IS JUST THE SAME AS  unemp  MINUS TREND MINUS SEASONAL

   plot(m1$data[,"seasonal"], ylab="Seasonal Unemp", main=
        "Seasonal Component of Unemployment, from `seasonal' package")
                            ### saved as   UnempPic5.pdf

   m2 = seas(AirPassengers)
   plot(m2)
   plot(m2$data[,"seasonal"], ylab="Seasonal Air Passengers")

### In both of these pictures, we want to understand hoiw to "predict" 
##      the seasonal part of a Seasonal ARIMA model

#-----------------------------------------------------------------------
# First the identification of a differenced and seasonally differenced
#   Seasonal Arima Model

 dfSunmp = diff(unemp) - lag(diff(unemp),12)    ## of length 323-13 = 310
 par(mfrow=c(2,1))
 acf(dfSunmp)
 pacf(dfSunmp)

 fit.unemp0 = arima(unemp, order=c(1,1,1), seasonal=list(order=c(0,1,1), period=12))
 fit.dfSunmp = arima(unemp-m1$data[,"trend"], order=c(1,0,1), include.mean=F,
         seasonal=list(order=c(0,0,1), period=12))

> fit.unemp0
Call:
arima(x = unemp, order = c(1, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12))

Coefficients:
         ar1      ma1     sma1
      0.9436  -0.8254  -0.8507
s.e.  0.0293   0.0430   0.0388

sigma^2 estimated as 62112:  log likelihood = -2157.98,  aic = 4323.96

### same model fit as that reported from "seasonal"

#------------------------------------------------------------------

> plot(fit.unemp0$residuals,m1$data[,"irregular"])
        ## generally increase together, but not a very close linear relationship

> unemp.fit = unemp-fit.unemp0$residuals
  plot(unemp)
  lines(unemp.fit, col="blue")

> par(mfrow=c(3,1))
  plot(unemp-m1$data[,"trend"], ylab="Detrended unemp", main=
      "unemp  minus  trend from package 'seasonal'")
  plot(m1$data[,"seasonal"], ylab="Seasonal unemp", main=
      "Seasonal component of  unemp  from package 'seasonal'")
  plot(unemp.fit-m1$data[,"trend"], ylab="Detrended", main=
      "1-Step-Ahead fitted  unemp  minus trend from `seasonal' package")
                     ### saved as   UnempPic6.pdf

  ###  SO THE PROCESS OF DECOMPOSING THE DETRENDED SERIES INTO SEASONAL AND IRREGULAR 
INVOLVES A LITTLE BIT OF SMOOTHING WITH THE EFFECT OF MAKING THE WITHIN-YEAR WAVE-FORMS 
SMOOTHER MORE SIMILAR TO ONE ANOTHER.

 ### HOW TO SEPARATE AND PROPAGATE THE WITHIN-YEAR WAVE-FORM ??
 ### NOTE THAT THE SEASONAL COMPONENT AND DETRENDED SERIES BOTH HAVE LARGEST PEAKS
 ###    SLOWLY SHRINKING YEAR BY YEAR WHILE THE SECONDARY PEAKS ARE GROWING,
 ###    BUT THE FITTED STEP-AHEAD SERIES DOES NOT SHOW THIS !!