set.seed(154)
w <- rnorm(100,0,1)
x <- filter(w, filter=c(0,-.9), method="recursive")
v <- filter(x, rep(1/4, 4), sides = 1)
plot.ts(x, main = "autoregression")
lines(v, col=10, lty=2)
cs <- cos(2*pi*(1:100)/4)
v2 <- filter(cs, rep(1/4, 4), sides = 1)
plot.ts(cs, main = expression(cos(2*pi*t/4)))
lines(v2, col=10, lty=2)
cs <- cos(2*pi*(1:100)/4); w = rnorm(100,0,1)
v3 <- filter(cs+w, rep(1/4, 4), sides = 1)
plot.ts(cs+w, main = expression(cos(2*pi*t/4)+N(0,1)))
lines(v3, col=10, lty=2)
In 1.3(a), we cannot find an obvious trend from the plot after taking moving averages while the series becomes mean-reverting along x=0.
Taking moving averages of a cosine function in 1.3(b) results in a horizontal line (x=0) because there’s no noise in the function.
We repeat the process in 1.3(b) with added N(0,1) noise in 1.3(c). Taking moving averages reduces the noise in original time series data, and the series becomes mean-reverting as in 1.3(a).
\(x_t=β_1+β_2t+w_t\)
According to \(E(w_t)=0\), we’ll have the mean function \(E(x_t)=E(β_1+β_2t+w_t)=β_1+β_2t\).
The mean function varies along time, which isn’t constant and independent of time. As a result, \(x_t\) is not stationary.
Firstly, since \(E(w_t)=0\), we’ll have the mean function of \(y_t\):
\(E(y_t)=E(x_t)-E(x_{t-1})=E(β_1+β_2t+w_t)-E(β_1+β_2(t-1)+w_{t-1})=β_2\)
The mean of \(y_t\) is a constant.
Secondly, the auto-covariance \(\gamma(s,t)=Cov(y_s,y_t)\)
When s = t the autocovariance of \(y_t\) is
\(Cov(y_t,y_t)\)
\(=Cov(x_t-x_{t-1},x_t-x_{t-1})\)
\(=Cov(β_2+w_t-w_{t-1},β_2+w_t-w_{t-1})\)
\(=Cov(w_t,w_t)+Cov(w_{t-1},w_{t-1})\)
\(=2σ_{w}^2\)
When s = t+1 \(Cov(x_{t+1}-x_t,x_t-x_{t-1})\)
\(=Cov((β_1+β_2(t+1)+w_{t+1})-(β_1+β_2t+w_t),(β_1+β_2t+w_t)-(β_1+β_2(t-1)+w_{t-1}))\)
\(=Cov(β_2+w_{t+1}-w_t,β_2+w_t-w_{t-1})\)
\(=-Cov(w_t,w_t)\)
\(=-σ_{w}^2\)
When s = t-1
\(Cov(y_{t-1},y_t)\)
\(=Cov(x_{t-1}-x_{t-2},x_t-x_{t-1})\)
\(=Cov((β_1+β_2(t-1)+w_{t-1})-(β_1+β_2(t-2)+w_{t-2}),β_2+w_t-w_{t-1})\)
\(=Cov(β_2+w_{t-1}-w_{t-2},β_2+w_t-w_{t-1})\)
\(=-σ_{w}^2\)
As a result,\(\gamma_{y}(r,s)\):
\((i)\) \(2σ_{w}^2, if h=0\)
\((ii)\) \(−σ_{w}^2, if h=1\)
\((iii)\) \(0, otherwise\)
\(x_t=β_1+β_2t+w_t\),
\(V_t=\frac{1}{2q+1} \sum_{j=-q}^q x_{t-j}\)
\(E(V_t)=E\left(\frac{1}{2q+1}\sum_{j=-q}^q X_{t-j}\right)\)
\(=\frac{1}{2q+1}E\left(x_{t+q}+x_{t+q-1}+...+x_{t-q}\right)\)
\(=\frac{1}{2q+1}E\left((β_1+β_2(t+q)+w_{t+q})+...+(β_1+β_2(t-q)+w_{t-q})\right)\)
\(=\frac{1}{2q+1}E\left(β_1(2q+1)+β_2t(2q+1)+\sum_{j=-q}^q w_{t-j} \right)\)
Since \(E(w_t)=0\), we will have \(E(V_t)=β_1+β_2t\)
\(x_t=w_{t-1}+2w_t+w_{t+1}\)
Autocovariance \(\gamma(r,s)=Cov(x_s,x_t)\)
When s = t
\(Cov(x_t,x_t)\)
\(=Cov(w_{t-1}+2w_t+w_{t+1},w_{t-1}+2w_t+w_{t+1})\)
\(=Cov(w_{t-1},w_{t-1})+4Cov(w_t,w_t)+Cov(w_{t+1},w_{t+1})\)
\(=6σ_{w}^2\)
When s = t+1
\(Cov(x_{t+1},x_t)\)
\(=Cov(w_t+2w_{t+1}+w_{t+2},w_{t-1}+2w_t+w_{t+1})\)
\(=Cov(w_t,2w_t)+2Cov(w_{t+1},w_{t+1})\)
\(=4σ_{w}^2\)
When s = t+2
\(Cov(x_{t+2},x_t)\)
\(=Cov(w_{t+1}+2w_{t+2}+w_{t+3},w_{t-1}+2w_t+w_{t+1})\)
\(=Cov(w_{t+1},w_{t+1})\)
\(=σ_{w}^2\)
As a result,
\(\gamma_x(r,s)=\)
\((i)\) \(6σ_{w}^2, if s=t\)
\((ii)\) \(4σ_{w}^2, if \mid s-t\mid=1\)
\((iii)\) \(σ_{w}^2, if \mid s-t\mid=2\)
\((iiii)\) \(0, otherwise\)
The ACF formula:
\(\rho(s,t)=\frac{\gamma(s,t)}{\sqrt{\gamma(s,s)\gamma(t,t)}}\)
When s=t, \(\rho_x(t,t) = \frac{\gamma_x(t,t)}{\sqrt{\gamma_x(t,t) \gamma_x(t,t)}} = \frac{6σ_{w}^2}{\sqrt{(6σ_{w}^2)(6σ_{w}^2)}}=1\)
When s = t+1, \(\rho_x(t+1,t) = \frac{\gamma(t+1,t)}{\sqrt{\gamma_x(t+1,t+1)\gamma_x(t,t)}} = \frac{4σ_{w}^2}{\sqrt{(6σ_{w}^2)(6σ_{w}^2)}} = \frac{2}{3}\)
When s = t+2, \(\rho_x(t+2,t) = \frac{\gamma(t+2,t+2)}{\sqrt{\gamma_x(t+2,t+2)\gamma_x(t,t)}} = \frac{σ_{w}^2}{\sqrt{(6σ_{w}^2)(6σ_{w}^2)}} = \frac{1}{6}\)
\(\rho_x(s,t)=\)
\(1\) , if \(s=t\)
\(\frac{2}{3}\), if \(\mid s-t\mid=1\)
\(\frac{1}{6}\), if \(\mid s-t\mid=2\)
\(0\), otherwise
The ACF of \(x_t\)
w1 = rnorm(10000,0,1)
x1 = filter(w1, sides = 2, filter=c(1,2,1))
acf(x1[!is.na(x1)], lag.max = 5, main = " ")
\(x_t=w_t\), \(w_t satisfies N(0,σ_{w}^2)\)
\(y_t=w_t-\theta w_{t-1}+u_t\), \(u_t satisfies N(0,σ_{u}^2)\)
When s = t,
\(\gamma_y(t,t)\)
\(=Cov(w_t-\theta w_{t-1}+u_t,w_t-\theta w_{t-1}+u_t)\)
\(=Cov(w_t,w_t)+\theta^2 Cov(w_{t-1},w_{t-1})+Cov(u_t,u_t)\)
\(=\sigma^2_w+\theta^2 \sigma^2_w+\sigma^2_u\)
\(=(\theta^2+1)\sigma^2_w+\sigma^2_u\)
As a result,
\(\rho_y(s,t) = \rho_y(t,t) = \frac{\gamma_y(t,t)}{\sqrt{\gamma_y(t,t)\gamma_y(t,t)}}\)
\(\frac{(\theta^2+1)\sigma^2_w+\sigma^2_u}{\sqrt{((\theta^2+1)\sigma^2_w+\sigma^2_u)((\theta^2+1)\sigma^2_w+\sigma^2_u)}}\)
\(=1\)
When s = t+1
\(\gamma_y(t+1,t)\)
\(=Cov(w_{t+1}-\theta w_t+u_{t+1},w_t-\theta w_{t-1}+u_t)\)
\(=-\theta Cov(w_t,w_t)\)
\(=-\theta \sigma^2_w\)
\(\gamma_y(t+1,t+1)\)
\(=Cov(w_{t+1},w_{t+1})+\theta^2 Cov(w_t,w_t)+Cov(u_{t+1},u_{t+1})\)
\(=(\theta^2+1)\sigma^2_w+\sigma^2_u\)
As a result,
\(\rho_y(s,t)=\rho_y(t+1,t)=\frac{\gamma_y(t+1,t)}{\sqrt{\gamma_y(t+1,t+1)\gamma_y(t,t)}}\)
\(=\frac{-\theta \sigma^2_w}{\sqrt{((\theta^2+1)\sigma^2_w+\sigma^2_u)^2}}\)
\(=\frac{-\theta \sigma^2_w}{(\theta^2+1)\sigma^2_w+\sigma^2_u}\)
The ACF of \(y_t\) is
\(1\) , if \(s = t\)
\(\frac{-\theta \sigma^2_w}{(\theta^2+1)\sigma^2_w+\sigma^2_u}\), if \(\mid s-t\mid = 1\)
\(0\) , \(otherwise\)
The CCF formula \(= \rho_xy(h) = \rho_xy(s,t) = \frac{\gamma_xy(s,t)}{\sqrt{\gamma_x(s,s)\gamma_y(t,t)}}\)
When s = t
\(\gamma_x(t,t)\)
\(=Cov(x_t,x_t)\)
\(=Cov(w_t,w_t)\)
\(=\sigma^2_w\)
\(\gamma_y(t,t)\)
\(=Cov(y_t,y_t)\)
\(=(\theta^2+1)\sigma^2_w+\sigma^2_u\)
\(\gamma_xy(x_t,y_t)\)
\(=Cov(w_t,w_t-\theta w_{t-1}+u_t)\)
\(=Cov(w_t,w_t)\)
\(=\sigma^2_w\)
\(\rho_xy(t,t)=\frac{Cov(x_t,y_t)}{\sqrt{\gamma_x(t,t)\gamma_y(t,t)}}\)
\(=\frac{\sigma^2_w}{\sqrt{(\theta^2+1)\sigma^4_w+\sigma^2_w\sigma^2_u}}\)
When s = t+1
\(\rho_xy(t+1,t)\)
\(=\frac{Cov(x_{t+1},y_t)}{\sqrt{\gamma_x(t+1,t+1)\gamma_y(t,t)}}\)
Since \(Cov(x_{t+1},y_t) = Cov(w_{t+1},w_t-\theta w_{t-1}+u_t) = 0\)
\(\rho_xy(t+1,t) = 0\)
When s = t-1
\(Cov(x_{t-1},y_t) = Cov(w_{t-1},w_t-\theta w_{t-1}+u_t) = -\theta \sigma^2_w\)
\(\rho_xy(t-1,t)\)
\(=\frac{-theta \sigma^2_w}{\sqrt{(\theta^2+1)\sigma^4_w+\sigma^2_w \sigma^2_u}}\)
The CCF of the two series \(x_t\) and \(y_t\) is
\(\frac{\sigma^2_w}{\sqrt{(\theta^2+1)\sigma^4_w+\sigma^2_w\sigma^2_u}}\) , if \(s = t\)
\(0\), if \(s=t+1\)
\(\frac{-\theta \sigma^2_w}{\sqrt{(\theta^2+1)\sigma^4_w+\sigma^2_w \sigma^2_u}}\) , if \(s=t-1\)
\(0\) , \(otherwise\)
require(astsa)
## Loading required package: astsa
trend = time(jj) - 1970
Q = factor(cycle(jj))
reg = lm(log(jj)~0 + trend + Q, na.action=NULL)
summary(reg)
##
## Call:
## lm(formula = log(jj) ~ 0 + trend + Q, na.action = NULL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.29318 -0.09062 -0.01180 0.08460 0.27644
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## trend 0.167172 0.002259 74.00 <2e-16 ***
## Q1 1.052793 0.027359 38.48 <2e-16 ***
## Q2 1.080916 0.027365 39.50 <2e-16 ***
## Q3 1.151024 0.027383 42.03 <2e-16 ***
## Q4 0.882266 0.027412 32.19 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1254 on 79 degrees of freedom
## Multiple R-squared: 0.9935, Adjusted R-squared: 0.9931
## F-statistic: 2407 on 5 and 79 DF, p-value: < 2.2e-16
trend = time(jj) - 1970
Q = factor(cycle(jj))
reg = lm(log(jj)~1 + trend + Q, na.action=NULL)
summary(reg)
##
## Call:
## lm(formula = log(jj) ~ 1 + trend + Q, na.action = NULL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.29318 -0.09062 -0.01180 0.08460 0.27644
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.052793 0.027359 38.480 < 2e-16 ***
## trend 0.167172 0.002259 73.999 < 2e-16 ***
## Q2 0.028123 0.038696 0.727 0.4695
## Q3 0.098231 0.038708 2.538 0.0131 *
## Q4 -0.170527 0.038729 -4.403 3.31e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1254 on 79 degrees of freedom
## Multiple R-squared: 0.9859, Adjusted R-squared: 0.9852
## F-statistic: 1379 on 4 and 79 DF, p-value: < 2.2e-16
fit1 = lm(log(jj)~0 + trend + Q, na.action=NULL)
ts.plot(cbind(log(jj),fit1$fitted),col=1:2)
fit2 = lm(log(jj)~1 + trend + Q, na.action=NULL)
ts.plot(cbind(log(jj),fit2$fitted),col=3:4)
ts.plot(fit1$resid)
ts.plot(fit2$resid)
par(mfrow=c(3,1))
plot(cmort, main="Cardiovascular Mortality", xlab="", ylab="")
plot(tempr, main="Temperature", xlab="", ylab="")
plot(part, main="Particulates", xlab="", ylab="")
pairs(cbind(Mortality=cmort, Temperature=tempr, Particulates=part))
temp = tempr-mean(tempr)
temp2 = temp^2
trend = time(cmort)
partL4=lag(part,-4)
newdata = ts.intersect(cmort, trend, temp, temp2, part, partL4, dframe=TRUE)
summary(lm(cmort~ trend + temp + temp2 + part + partL4, data = newdata, na.action = NULL))
##
## Call:
## lm(formula = cmort ~ trend + temp + temp2 + part + partL4, data = newdata,
## na.action = NULL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.228 -4.314 -0.614 3.713 27.800
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.808e+03 1.989e+02 14.123 < 2e-16 ***
## trend -1.385e+00 1.006e-01 -13.765 < 2e-16 ***
## temp -4.058e-01 3.528e-02 -11.503 < 2e-16 ***
## temp2 2.155e-02 2.803e-03 7.688 8.02e-14 ***
## part 2.029e-01 2.266e-02 8.954 < 2e-16 ***
## partL4 1.030e-01 2.485e-02 4.147 3.96e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.287 on 498 degrees of freedom
## Multiple R-squared: 0.608, Adjusted R-squared: 0.6041
## F-statistic: 154.5 on 5 and 498 DF, p-value: < 2.2e-16
par(mfrow=c(2,2))
plot(cmort, main="Cardiovascular Mortality", xlab="", ylab="")
plot(tempr, main="Temperature", xlab="", ylab="")
plot(part, main="Particulates", xlab="", ylab="")
plot(partL4, main="Lag Particulates", xlab="", ylab="")
pairs(cbind(Mortality=cmort, Temperature=tempr, Particulates=part, LagParticulates=partL4))
cor(cbind(cmort=newdata[,1],temp=newdata[,3],part=newdata[,5],partL4=newdata[,6]))
## cmort temp part partL4
## cmort 1.0000000 -0.4369648 0.4422896 0.5209993
## temp -0.4369648 1.0000000 -0.0148241 -0.3990848
## part 0.4422896 -0.0148241 1.0000000 0.5340505
## partL4 0.5209993 -0.3990848 0.5340505 1.0000000
temp = tempr-mean(tempr)
temp2 = temp^2
trend = time(cmort)
fit2.18 = lm(cmort~ trend, na.action = NULL)
fit2.19 = lm(cmort~ trend + temp, na.action = NULL)
fit2.20 = lm(cmort~ trend + temp + temp2, na.action = NULL)
fit2.21 = lm(cmort~ trend + temp + temp2 + part, na.action = NULL)
summary(fit2.18)
##
## Call:
## lm(formula = cmort ~ trend, na.action = NULL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16.445 -6.670 -1.366 5.505 40.107
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3297.6062 276.3132 11.93 <2e-16 ***
## trend -1.6249 0.1399 -11.61 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.893 on 506 degrees of freedom
## Multiple R-squared: 0.2104, Adjusted R-squared: 0.2089
## F-statistic: 134.9 on 1 and 506 DF, p-value: < 2.2e-16
summary(aov(fit2.18))
## Df Sum Sq Mean Sq F value Pr(>F)
## trend 1 10667 10667 134.9 <2e-16 ***
## Residuals 506 40020 79
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(aov(lm(cmort~cbind(trend))))
## Df Sum Sq Mean Sq F value Pr(>F)
## cbind(trend) 1 10667 10667 134.9 <2e-16 ***
## Residuals 506 40020 79
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(fit2.19)
##
## Call:
## lm(formula = cmort ~ trend + temp, na.action = NULL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16.846 -5.330 -1.207 4.701 33.306
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3125.75988 245.48233 12.73 <2e-16 ***
## trend -1.53785 0.12430 -12.37 <2e-16 ***
## temp -0.45792 0.03893 -11.76 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.887 on 505 degrees of freedom
## Multiple R-squared: 0.3802, Adjusted R-squared: 0.3778
## F-statistic: 154.9 on 2 and 505 DF, p-value: < 2.2e-16
summary(aov(fit2.19))
## Df Sum Sq Mean Sq F value Pr(>F)
## trend 1 10667 10667 171.5 <2e-16 ***
## temp 1 8607 8607 138.4 <2e-16 ***
## Residuals 505 31413 62
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(aov(lm(cmort~cbind(trend, temp))))
## Df Sum Sq Mean Sq F value Pr(>F)
## cbind(trend, temp) 2 19273 9637 154.9 <2e-16 ***
## Residuals 505 31413 62
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(fit2.20)
##
## Call:
## lm(formula = cmort ~ trend + temp + temp2, na.action = NULL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.464 -4.858 -0.945 4.511 34.939
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.038e+03 2.322e+02 13.083 < 2e-16 ***
## trend -1.494e+00 1.176e-01 -12.710 < 2e-16 ***
## temp -4.808e-01 3.689e-02 -13.031 < 2e-16 ***
## temp2 2.583e-02 3.287e-03 7.858 2.38e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.452 on 504 degrees of freedom
## Multiple R-squared: 0.4479, Adjusted R-squared: 0.4446
## F-statistic: 136.3 on 3 and 504 DF, p-value: < 2.2e-16
summary(aov(fit2.20))
## Df Sum Sq Mean Sq F value Pr(>F)
## trend 1 10667 10667 192.11 < 2e-16 ***
## temp 1 8607 8607 155.00 < 2e-16 ***
## temp2 1 3429 3429 61.75 2.38e-14 ***
## Residuals 504 27985 56
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(aov(lm(cmort~cbind(trend, temp, temp2))))
## Df Sum Sq Mean Sq F value Pr(>F)
## cbind(trend, temp, temp2) 3 22702 7567 136.3 <2e-16 ***
## Residuals 504 27985 56
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(fit2.21)
##
## Call:
## lm(formula = cmort ~ trend + temp + temp2 + part, na.action = NULL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.0760 -4.2153 -0.4878 3.7435 29.2448
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.831e+03 1.996e+02 14.19 < 2e-16 ***
## trend -1.396e+00 1.010e-01 -13.82 < 2e-16 ***
## temp -4.725e-01 3.162e-02 -14.94 < 2e-16 ***
## temp2 2.259e-02 2.827e-03 7.99 9.26e-15 ***
## part 2.554e-01 1.886e-02 13.54 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.385 on 503 degrees of freedom
## Multiple R-squared: 0.5954, Adjusted R-squared: 0.5922
## F-statistic: 185 on 4 and 503 DF, p-value: < 2.2e-16
summary(aov(fit2.21))
## Df Sum Sq Mean Sq F value Pr(>F)
## trend 1 10667 10667 261.62 <2e-16 ***
## temp 1 8607 8607 211.09 <2e-16 ***
## temp2 1 3429 3429 84.09 <2e-16 ***
## part 1 7476 7476 183.36 <2e-16 ***
## Residuals 503 20508 41
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(aov(lm(cmort~cbind(trend, temp, temp2, part))))
## Df Sum Sq Mean Sq F value Pr(>F)
## cbind(trend, temp, temp2, part) 4 30178 7545 185 <2e-16 ***
## Residuals 503 20508 41
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
num = length(cmort)
AIC(fit2.18)/num - log(2*pi)
## [1] 5.37846
BIC(fit2.18)/num - log(2*pi)
## [1] 5.403443
AIC(fit2.19)/num - log(2*pi)
## [1] 5.14025
BIC(fit2.19)/num - log(2*pi)
## [1] 5.173561
AIC(fit2.20)/num - log(2*pi)
## [1] 5.028611
BIC(fit2.20)/num - log(2*pi)
## [1] 5.070249
AIC(fit2.21)/num - log(2*pi)
## [1] 4.721732
BIC(fit2.21)/num - log(2*pi)
## [1] 4.771699