HW2
ts(1:10, frequency = 4, start = c(1959, 2)) # 2nd Quarter of 1959
## Qtr1 Qtr2 Qtr3 Qtr4
## 1959 1 2 3
## 1960 4 5 6 7
## 1961 8 9 10
makeDF = function(x) {
df = data.frame(
y = log10(x),
t = time(x),
s = as.factor(cycle(x))
)
}
seaLS1 = function(x) {
df = makeDF(x)
lm(y ~ t - 1 + s, data = df)
}
seaLS2 = function(x) {
df = makeDF(x)
lm(y ~ t + s, data = df, contrasts = list(s = contr.helmert))
}
mod1 = seaLS1(JohnsonJohnson)
mod2 = seaLS2(JohnsonJohnson)
model.matrix(mod1)
## t s1 s2 s3 s4
## 1 1960.00 1 0 0 0
## 2 1960.25 0 1 0 0
## 3 1960.50 0 0 1 0
## 4 1960.75 0 0 0 1
## 5 1961.00 1 0 0 0
## 6 1961.25 0 1 0 0
## 7 1961.50 0 0 1 0
## 8 1961.75 0 0 0 1
## 9 1962.00 1 0 0 0
## 10 1962.25 0 1 0 0
## 11 1962.50 0 0 1 0
## 12 1962.75 0 0 0 1
## 13 1963.00 1 0 0 0
## 14 1963.25 0 1 0 0
## 15 1963.50 0 0 1 0
## 16 1963.75 0 0 0 1
## 17 1964.00 1 0 0 0
## 18 1964.25 0 1 0 0
## 19 1964.50 0 0 1 0
## 20 1964.75 0 0 0 1
## 21 1965.00 1 0 0 0
## 22 1965.25 0 1 0 0
## 23 1965.50 0 0 1 0
## 24 1965.75 0 0 0 1
## 25 1966.00 1 0 0 0
## 26 1966.25 0 1 0 0
## 27 1966.50 0 0 1 0
## 28 1966.75 0 0 0 1
## 29 1967.00 1 0 0 0
## 30 1967.25 0 1 0 0
## 31 1967.50 0 0 1 0
## 32 1967.75 0 0 0 1
## 33 1968.00 1 0 0 0
## 34 1968.25 0 1 0 0
## 35 1968.50 0 0 1 0
## 36 1968.75 0 0 0 1
## 37 1969.00 1 0 0 0
## 38 1969.25 0 1 0 0
## 39 1969.50 0 0 1 0
## 40 1969.75 0 0 0 1
## 41 1970.00 1 0 0 0
## 42 1970.25 0 1 0 0
## 43 1970.50 0 0 1 0
## 44 1970.75 0 0 0 1
## 45 1971.00 1 0 0 0
## 46 1971.25 0 1 0 0
## 47 1971.50 0 0 1 0
## 48 1971.75 0 0 0 1
## 49 1972.00 1 0 0 0
## 50 1972.25 0 1 0 0
## 51 1972.50 0 0 1 0
## 52 1972.75 0 0 0 1
## 53 1973.00 1 0 0 0
## 54 1973.25 0 1 0 0
## 55 1973.50 0 0 1 0
## 56 1973.75 0 0 0 1
## 57 1974.00 1 0 0 0
## 58 1974.25 0 1 0 0
## 59 1974.50 0 0 1 0
## 60 1974.75 0 0 0 1
## 61 1975.00 1 0 0 0
## 62 1975.25 0 1 0 0
## 63 1975.50 0 0 1 0
## 64 1975.75 0 0 0 1
## 65 1976.00 1 0 0 0
## 66 1976.25 0 1 0 0
## 67 1976.50 0 0 1 0
## 68 1976.75 0 0 0 1
## 69 1977.00 1 0 0 0
## 70 1977.25 0 1 0 0
## 71 1977.50 0 0 1 0
## 72 1977.75 0 0 0 1
## 73 1978.00 1 0 0 0
## 74 1978.25 0 1 0 0
## 75 1978.50 0 0 1 0
## 76 1978.75 0 0 0 1
## 77 1979.00 1 0 0 0
## 78 1979.25 0 1 0 0
## 79 1979.50 0 0 1 0
## 80 1979.75 0 0 0 1
## 81 1980.00 1 0 0 0
## 82 1980.25 0 1 0 0
## 83 1980.50 0 0 1 0
## 84 1980.75 0 0 0 1
## attr(,"assign")
## [1] 1 2 2 2 2
## attr(,"contrasts")
## attr(,"contrasts")$s
## [1] "contr.treatment"
model.matrix(mod2)
## (Intercept) t s1 s2 s3
## 1 1 1960.00 -1 -1 -1
## 2 1 1960.25 1 -1 -1
## 3 1 1960.50 0 2 -1
## 4 1 1960.75 0 0 3
## 5 1 1961.00 -1 -1 -1
## 6 1 1961.25 1 -1 -1
## 7 1 1961.50 0 2 -1
## 8 1 1961.75 0 0 3
## 9 1 1962.00 -1 -1 -1
## 10 1 1962.25 1 -1 -1
## 11 1 1962.50 0 2 -1
## 12 1 1962.75 0 0 3
## 13 1 1963.00 -1 -1 -1
## 14 1 1963.25 1 -1 -1
## 15 1 1963.50 0 2 -1
## 16 1 1963.75 0 0 3
## 17 1 1964.00 -1 -1 -1
## 18 1 1964.25 1 -1 -1
## 19 1 1964.50 0 2 -1
## 20 1 1964.75 0 0 3
## 21 1 1965.00 -1 -1 -1
## 22 1 1965.25 1 -1 -1
## 23 1 1965.50 0 2 -1
## 24 1 1965.75 0 0 3
## 25 1 1966.00 -1 -1 -1
## 26 1 1966.25 1 -1 -1
## 27 1 1966.50 0 2 -1
## 28 1 1966.75 0 0 3
## 29 1 1967.00 -1 -1 -1
## 30 1 1967.25 1 -1 -1
## 31 1 1967.50 0 2 -1
## 32 1 1967.75 0 0 3
## 33 1 1968.00 -1 -1 -1
## 34 1 1968.25 1 -1 -1
## 35 1 1968.50 0 2 -1
## 36 1 1968.75 0 0 3
## 37 1 1969.00 -1 -1 -1
## 38 1 1969.25 1 -1 -1
## 39 1 1969.50 0 2 -1
## 40 1 1969.75 0 0 3
## 41 1 1970.00 -1 -1 -1
## 42 1 1970.25 1 -1 -1
## 43 1 1970.50 0 2 -1
## 44 1 1970.75 0 0 3
## 45 1 1971.00 -1 -1 -1
## 46 1 1971.25 1 -1 -1
## 47 1 1971.50 0 2 -1
## 48 1 1971.75 0 0 3
## 49 1 1972.00 -1 -1 -1
## 50 1 1972.25 1 -1 -1
## 51 1 1972.50 0 2 -1
## 52 1 1972.75 0 0 3
## 53 1 1973.00 -1 -1 -1
## 54 1 1973.25 1 -1 -1
## 55 1 1973.50 0 2 -1
## 56 1 1973.75 0 0 3
## 57 1 1974.00 -1 -1 -1
## 58 1 1974.25 1 -1 -1
## 59 1 1974.50 0 2 -1
## 60 1 1974.75 0 0 3
## 61 1 1975.00 -1 -1 -1
## 62 1 1975.25 1 -1 -1
## 63 1 1975.50 0 2 -1
## 64 1 1975.75 0 0 3
## 65 1 1976.00 -1 -1 -1
## 66 1 1976.25 1 -1 -1
## 67 1 1976.50 0 2 -1
## 68 1 1976.75 0 0 3
## 69 1 1977.00 -1 -1 -1
## 70 1 1977.25 1 -1 -1
## 71 1 1977.50 0 2 -1
## 72 1 1977.75 0 0 3
## 73 1 1978.00 -1 -1 -1
## 74 1 1978.25 1 -1 -1
## 75 1 1978.50 0 2 -1
## 76 1 1978.75 0 0 3
## 77 1 1979.00 -1 -1 -1
## 78 1 1979.25 1 -1 -1
## 79 1 1979.50 0 2 -1
## 80 1 1979.75 0 0 3
## 81 1 1980.00 -1 -1 -1
## 82 1 1980.25 1 -1 -1
## 83 1 1980.50 0 2 -1
## 84 1 1980.75 0 0 3
## attr(,"assign")
## [1] 0 1 2 2 2
## attr(,"contrasts")
## attr(,"contrasts")$s
## [,1] [,2] [,3]
## 1 -1 -1 -1
## 2 1 -1 -1
## 3 0 2 -1
## 4 0 0 3
summary(mod1)
##
## Call:
## lm(formula = y ~ t - 1 + s, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.127327 -0.039356 -0.005123 0.036740 0.120056
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## t 7.260e-02 9.811e-04 74.00 <2e-16 ***
## s1 -1.426e+02 1.933e+00 -73.76 <2e-16 ***
## s2 -1.426e+02 1.933e+00 -73.75 <2e-16 ***
## s3 -1.425e+02 1.933e+00 -73.72 <2e-16 ***
## s4 -1.426e+02 1.934e+00 -73.77 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05445 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
summary(mod2)
##
## Call:
## lm(formula = y ~ t + s, data = df, contrasts = list(s = contr.helmert))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.127327 -0.039356 -0.005123 0.036740 0.120056
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.426e+02 1.933e+00 -73.751 < 2e-16 ***
## t 7.260e-02 9.811e-04 73.999 < 2e-16 ***
## s1 6.107e-03 8.403e-03 0.727 0.4695
## s2 1.218e-02 4.852e-03 2.511 0.0141 *
## s3 -2.309e-02 3.432e-03 -6.727 2.48e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05445 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
library(car)
## Loading required package: carData
linearHypothesis(mod1, "s3=s4")
## Linear hypothesis test
##
## Hypothesis:
## s3 - s4 = 0
##
## Model 1: restricted model
## Model 2: y ~ t - 1 + s
##
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 80 0.37724
## 2 79 0.23422 1 0.14302 48.238 9.503e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(log10(JohnsonJohnson))
x = JohnsonJohnson
window(x) = predict(mod1)
lines(x, col = "red")
