Question 1)
Example 1.2
Part A
# Defining the sample size
n <- 200
# Generating the signal st for both parts of the interval, from 0 to 100 observations we would get 0
# for 101 to 100 observations we would get a value according to the function
st <- c(rep(0, 100), 10 * exp(-(1:100 - 100) / 20) * cos(2 * pi * (1:100) / 4))
# Generate thing Gaussian white noise wt with variance 1
wt <- rnorm(n, mean = 0, sd = 1)
# Generating the observed signal xt
xt <- st + wt
plot.ts(xt, main="Signal-plus-Noise Model Simulation for part A", ylab="Value", xlab="Time t")
B)
# Generating the signal st for both parts of the interval, from 0 to 100 observations we would get 0
# for 101 to 100 observations we would get a value according to the function
stb <- c(rep(0, 100), 10 * exp(-(1:100 - 100) / 200) * cos(2 * pi * (1:100) / 4))
# Generate thing Gaussian white noise wt with variance 1
wtb <- rnorm(n, mean = 0, sd = 1)
# Generating the observed signal xt
xtb <- stb + wtb
plot.ts(xtb, main="Signal-plus-Noise Model Simulation for part B", ylab="Value", xlab="Time t")
c
# Compare the signal modulators exp(-t/20) and exp(-t/200)
t_values <- 1:100
modulator_a <- exp(-t_values / 20)
modulator_b <- exp(-t_values / 200)
# Plotting the signal modulators with adjusted y-axis limits and increased line width
plot(t_values, modulator_a, type='l', col='blue', ylim=c(0, 1), xlab='Time t', ylab='Value',
main='Comparison of Signal Modulators', lwd=2)
lines(t_values, modulator_b, type='l', col='red', lwd=2)
# Adding a legend to the plot
legend("topright", inset=.05, legend=c('Part a- exp(-t/20)', 'Part b- exp(-t/200)'), col=c('blue', 'red'), lty=1, lwd=2, bty="n")
'
The simulated series (a) has a fast-decaying oscillation which might resemble the decaying
oscillations after the initial shock of an earthquake. Series (b) features a slower decay and
could potentially mimic the prolonged aftershocks of an earthquake. However, neither series
perfectly models the sharp spike followed by decay seen in the explosion series, as both (a)
and (b) include oscillatory components which are not present in the explosion data. The
explosion data would be better modeled by a single spike without periodic oscillations. The
white noise component in both series adds randomness, similar to the background noise observed
in the earthquake and explosion series.
'
## [1] "\nThe simulated series (a) has a fast-decaying oscillation which might resemble the decaying\noscillations after the initial shock of an earthquake. Series (b) features a slower decay and\ncould potentially mimic the prolonged aftershocks of an earthquake. However, neither series\nperfectly models the sharp spike followed by decay seen in the explosion series, as both (a)\nand (b) include oscillatory components which are not present in the explosion data. The\nexplosion data would be better modeled by a single spike without periodic oscillations. The\nwhite noise component in both series adds randomness, similar to the background noise observed\nin the earthquake and explosion series.\n"
Question 1.3
Part A
# Set seed for reproducibility
set.seed(154)
# Parameters
n <- 100 # Number of observations
sigma_w <- 1 # Standard deviation of w
# Generate white noise w
w <- rnorm(n, mean = 0, sd = sigma_w)
# Generate x using the autoregression model
x <- rep(0, n)
for (t in 3:n) {
x[t] <- -0.9 * x[t-2] + w[t]
}
# Apply moving average filter
v <- filter(x, rep(1/4, 4), sides = 1)
# Plotting
v <- na.omit(v) # Note: This might change the length of v
# Plotting
plot(x, type = 'l', col = 'blue', ylim = range(c(x, v, na.rm = TRUE)), xlab = "Time", ylab = "Value", main = "Autoregression and Moving Average Filter")
lines(v, col = 'red', lty = 2) # Adding the moving average as a dashed line
legend("topright", legend = c("x_t", "v_t (MA Filter)"), col = c("blue", "red"), lty = c(1, 2), cex = 0.8)
'
There is a significant autocorrelation at lag 2, as each value is directly dependent on the
value two steps before, modified by the effect of the noise.
When observing the plot, the smoothed series vt will typically show less variability than
xt, with extreme values and short-term fluctuations being averaged out.
This can be particularly useful for highlighting more gradual trends or patterns in the data
that might be obscured by noise in the original series.
The moving average thus serves as a tool for data analysis, helping to identify underlying
trends by reducing noise and making the series easier to analyze for patterns, trends, or
cycles.
'
## [1] "\nThere is a significant autocorrelation at lag 2, as each value is directly dependent on the\nvalue two steps before, modified by the effect of the noise.\nWhen observing the plot, the smoothed series vt will typically show less variability than \n xt, with extreme values and short-term fluctuations being averaged out.\n \nThis can be particularly useful for highlighting more gradual trends or patterns in the data\nthat might be obscured by noise in the original series.\nThe moving average thus serves as a tool for data analysis, helping to identify underlying\ntrends by reducing noise and making the series easier to analyze for patterns, trends, or\ncycles.\n"
t <- 1:n # Time index
x <- cos(2 * pi * t / 4) # Autoregression
# Apply moving average filter
v <- filter(x, rep(1/4, 4), sides = 1)
# Plotting
v <- na.omit(v) # Note: This might change the length of v
# Plotting
plot(x, type = 'l', col = 'blue', ylim = range(c(x, v, na.rm = TRUE)), xlab = "Time", ylab = "Value", main = "Autoregression and Moving Average Filter")
lines(v, col = 'red', lty = 2) # Adding the moving average as a dashed line
legend("topright", legend = c("x_t", "v_t (MA Filter)"), col = c("blue", "red"), lty = c(1, 2), cex = 0.8)
"
The original series xt demonstrates the predictable oscillatory behavior of a cosine function,
with regular peaks and troughs. The application of the moving average filter to create
vt modifies this behavior by smoothing the series, introducing a phase shift, and reducing
the amplitude of oscillations. These changes highlight the filter's effectiveness in reducing
noise and making underlying patterns more apparent
"
## [1] "\nThe original series xt demonstrates the predictable oscillatory behavior of a cosine function,\nwith regular peaks and troughs. The application of the moving average filter to create \nvt modifies this behavior by smoothing the series, introducing a phase shift, and reducing\nthe amplitude of oscillations. These changes highlight the filter's effectiveness in reducing\nnoise and making underlying patterns more apparent\n"
set.seed(154) # Ensure reproducibility
n <- 100 # Number of observations
t <- 1:n # Time index
w <- rnorm(n) # Generate N(0, 1) noise
x <- cos(2 * pi * t / 4) + w # Add noise to the cos wave
# Apply the moving average filter
v <- filter(x, rep(1/4, 4), sides = 1)
# Plotting
plot(x, type = 'l', col = 'blue', xlab = "Time", ylab = "Value", main = "Noisy Cosine Series and Moving Average Filter")
lines(v, col = 'red', lty = 2) # Add the moving average as a dashed line
legend("topright", legend = c("x_t", "v_t"), col = c("blue", "red"), lty = c(1, 2), cex = 0.8)
"
The addition of white noise wt to the cosine series xt increases its volatility and obscures
its clean periodic pattern, making the series more erratic. Applying the moving average filter
vt counteracts some of these effects by smoothing the series, reducing noise-induced
fluctuations, and making the underlying periodic pattern more apparent, albeit with reduced
amplitude and a slight phase shift. This demonstrates the utility of moving average filters in
extracting underlying patterns from noisy data, a common challenge in time series analysis.
"
## [1] "\nThe addition of white noise wt to the cosine series xt increases its volatility and obscures\nits clean periodic pattern, making the series more erratic. Applying the moving average filter \nvt counteracts some of these effects by smoothing the series, reducing noise-induced\nfluctuations, and making the underlying periodic pattern more apparent, albeit with reduced\namplitude and a slight phase shift. This demonstrates the utility of moving average filters in\nextracting underlying patterns from noisy data, a common challenge in time series analysis.\n"
"
Autoregressive Series: The filter smooths out the series, reducing the impact of noise and
making longer-term trends more visible, albeit at the cost of slightly obscuring the specific
autoregressive properties due to averaging out the lagged effects.
Pure Cosine Series: Given the series' inherent smoothness and regular oscillations, the moving
average filter's impact is subtle, slightly reducing the amplitude of the wave and introducing
a minor phase shift, but without significantly altering the series' fundamental periodic nature.
Cosine Series with Noise: The addition of white noise introduces significant volatility,
obscuring the underlying cosine pattern. The moving average filter is particularly effective
here, substantially reducing noise and revealing the cosine wave's periodicity with reduced
amplitude and a slight lag.
"
## [1] "\nAutoregressive Series: The filter smooths out the series, reducing the impact of noise and\nmaking longer-term trends more visible, albeit at the cost of slightly obscuring the specific\nautoregressive properties due to averaging out the lagged effects.\n\nPure Cosine Series: Given the series' inherent smoothness and regular oscillations, the moving\naverage filter's impact is subtle, slightly reducing the amplitude of the wave and introducing\na minor phase shift, but without significantly altering the series' fundamental periodic nature.\n\nCosine Series with Noise: The addition of white noise introduces significant volatility,\nobscuring the underlying cosine pattern. The moving average filter is particularly effective\nhere, substantially reducing noise and revealing the cosine wave's periodicity with reduced\namplitude and a slight lag.\n"
### Question 2)
## Part 1)
sales <- read.table("/Users/fahadmehfooz/Downloads/skisales.csv",header=TRUE,sep=",")
sales
## Quarter Sales PDI Season
## 1 Q1/64 37.0 109 1
## 2 Q2/64 33.5 115 0
## 3 Q3/64 30.8 113 0
## 4 Q4/64 37.9 116 1
## 5 Q1/65 37.4 118 1
## 6 Q2/65 31.6 120 0
## 7 Q3/65 34.0 122 0
## 8 Q4/65 38.1 124 1
## 9 Q1/66 40.0 126 1
## 10 Q2/66 35.0 128 0
## 11 Q3/66 34.9 130 0
## 12 Q4/66 40.2 132 1
## 13 Q1/67 41.9 133 1
## 14 Q2/67 34.7 135 0
## 15 Q3/67 38.8 138 0
## 16 Q4/67 43.7 140 1
## 17 Q1/68 44.2 143 1
## 18 Q2/68 40.4 147 0
## 19 Q3/68 38.4 148 0
## 20 Q4/68 45.4 151 1
## 21 Q1/69 44.9 153 1
## 22 Q2/69 41.6 156 0
## 23 Q3/69 44.0 160 0
## 24 Q4/69 48.1 163 1
## 25 Q1/70 49.7 166 1
## 26 Q2/70 43.9 171 0
## 27 Q3/70 41.6 174 0
## 28 Q4/70 51.0 175 1
## 29 Q1/71 52.0 180 1
## 30 Q2/71 46.2 184 0
## 31 Q3/71 47.1 187 0
## 32 Q4/71 52.7 189 1
## 33 Q1/72 52.2 191 1
## 34 Q2/72 47.0 193 0
## 35 Q3/72 47.8 194 0
## 36 Q4/72 52.8 196 1
## 37 Q1/73 54.1 199 1
## 38 Q2/73 49.5 201 0
## 39 Q3/73 49.5 202 0
## 40 Q4/73 54.3 204 1
## A linear pattern can be seen between Sales and PDI
plot(Sales~PDI, data=sales, pch=16)
sales_model = lm(Sales~PDI,data = sales)
abline(sales_model)
summary(sales_model)
##
## Call:
## lm(formula = Sales ~ PDI, data = sales)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.229 -2.686 0.554 2.728 4.454
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.39215 2.53942 4.88 1.93e-05 ***
## PDI 0.19791 0.01602 12.35 7.09e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.019 on 38 degrees of freedom
## Multiple R-squared: 0.8006, Adjusted R-squared: 0.7953
## F-statistic: 152.5 on 1 and 38 DF, p-value: 7.089e-15
"
The F stat value and its p value suggests indicates that the model
significantly predicts the outcome variable. As we reject the null hypothesis here.
"
## [1] "\nThe F stat value and its p value suggests indicates that the model\nsignificantly predicts the outcome variable. As we reject the null hypothesis here.\n"
Part 2)
# Plotting the ACF of the residuals
acf( residuals(sales_model), main="ACF of Regression Residuals", lag.max = 5)
'
Lag 1: In the first-order autocorrelation, the spike at lag 1 does not extend beyond the blue
dashed significance bounds, suggesting that there is no strong evidence of first-order
autocorrelation
'
## [1] "\nLag 1: In the first-order autocorrelation, the spike at lag 1 does not extend beyond the blue\ndashed significance bounds, suggesting that there is no strong evidence of first-order\nautocorrelation\n"
## Lets write a significance test for the same
library("car")
## Loading required package: carData
durbinWatsonTest(sales_model, max.lag=1, simulate=FALSE, alternative="two.sided")
## lag Autocorrelation D-W Statistic
## 1 -0.0008866674 1.968394
## Alternative hypothesis: rho0
"The The Autocorrelation is approximately -0.0008866674, which is very close to zero.
This suggests that there is almost no first-order autocorrelation detected in the residuals.
A value near to dw = 2 suggests the same thing which is no correlation.
This confirms ACF too which also suggested no correlation.
"
## [1] "The The Autocorrelation is approximately -0.0008866674, which is very close to zero.\nThis suggests that there is almost no first-order autocorrelation detected in the residuals.\nA value near to dw = 2 suggests the same thing which is no correlation.\nThis confirms ACF too which also suggested no correlation.\n"
Part 3)
'
Xt = beta0 + beta1 * Z1 + et
Xt is the dependent variable (e.g., ski sales in millions of dollars) at time t,
Zt is the independent variable (e.g., personal disposable income in billions of current dollars) at time t,
beta0 is the intercept term,
beta1 is the slope coefficient,
et is the error term at time t, which is assumed to be serially correlated.
The assumptions which are followed in linear regression are:
1. Linear relationship: There is a linear relationship between the dependent variable and the
independent variable.
2. Serial correlation in errors: Unlike the standard assumption of independence of errors in
ordinary least squares regression, here, the error terms et are assumed to be serially
correlated, typically modeled as:
rho is the autocorrelation coefficient for the first-order autocorrelation, and
ut is a white noise error term with mean zero and constant variance.
3. Homoscedasticity: The error terms have constant variance, var(ut) = sigma **2
4. No perfect multicollinearity: In the matrix of explanatory variables, no column is an
exact linear function of others.
5. Errors are normally distributed (for inference purposes): ut is normally distributed.
'
## [1] "\nXt = beta0 + beta1 * Z1 + et\n \nXt is the dependent variable (e.g., ski sales in millions of dollars) at time t,\nZt is the independent variable (e.g., personal disposable income in billions of current dollars) at time t,\nbeta0 is the intercept term,\nbeta1 is the slope coefficient, \net is the error term at time t, which is assumed to be serially correlated.\n\nThe assumptions which are followed in linear regression are:\n\n1. Linear relationship: There is a linear relationship between the dependent variable and the\nindependent variable.\n\n2. Serial correlation in errors: Unlike the standard assumption of independence of errors in\nordinary least squares regression, here, the error terms et are assumed to be serially\ncorrelated, typically modeled as:\n\nrho is the autocorrelation coefficient for the first-order autocorrelation, and \nut is a white noise error term with mean zero and constant variance.\n\n3. Homoscedasticity: The error terms have constant variance, var(ut) = sigma **2\n \n4. No perfect multicollinearity: In the matrix of explanatory variables, no column is an\nexact linear function of others.\n\n5. Errors are normally distributed (for inference purposes): ut is normally distributed.\n"
library(lmtest)
data_copy <- sales
# Step 1: Run the initial OLS regression on the copied data
sales_model1 <- lm(Sales ~ PDI, data = data_copy)
# Initialize variables for the iterative process
rho <- 0
prev_rho <- 1 # To store the rho from the previous iteration for comparison
converged <- FALSE
tolerance <- 1e-6 # Convergence tolerance
max_iterations <- 100 # Maximum allowed iterations
iteration <- 0
# Begin the iterative process
while (!converged && iteration < max_iterations) {
iteration <- iteration + 1
# Compute the residuals from the current model
residuals <- resid(sales_model1)
# Regress the lagged residuals on the current residuals to estimate rho
# Exclude the first observation to avoid NA values
lagged_residuals <- c(NA, residuals[-length(residuals)])
resid_model <- lm(residuals[-1] ~ lagged_residuals[-1] - 1) # No intercept
# Update rho if the change is significant
new_rho <- coef(resid_model)[[1]] # Extract the estimated coefficient
if (abs(new_rho - prev_rho) < tolerance) {
converged <- TRUE # Convergence achieved
} else {
prev_rho <- new_rho # Update previous rho for the next iteration
}
# Transform the variables using the estimated rho
data_copy$lagged_Sales <- c(NA, data_copy$Sales[-length(data_copy$Sales)])
data_copy$lagged_PDI <- c(NA, data_copy$PDI[-length(data_copy$PDI)])
data_transformed <- data_copy[-1,]
data_transformed$Sales_corrected <- data_transformed$Sales - new_rho * data_transformed$lagged_Sales
data_transformed$PDI_corrected <- data_transformed$PDI - new_rho * data_transformed$lagged_PDI
# Re-estimate the regression model with transformed variables
sales_model1 <- lm(Sales_corrected ~ PDI_corrected, data = data_transformed)
}
cat("Converged after", iteration, "iterations with rho =", new_rho, "\n")
## Converged after 8 iterations with rho = 0.01123638
print(summary(sales_model1))
##
## Call:
## lm(formula = Sales_corrected ~ PDI_corrected, data = data_transformed)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.2007 -2.7547 0.0238 2.8575 4.4526
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.48912 2.63764 4.356 0.000101 ***
## PDI_corrected 0.20234 0.01672 12.103 1.98e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.015 on 37 degrees of freedom
## Multiple R-squared: 0.7984, Adjusted R-squared: 0.7929
## F-statistic: 146.5 on 1 and 37 DF, p-value: 1.976e-14
"Converged after 8 iterations with rho = 0.01123638"
## [1] "Converged after 8 iterations with rho = 0.01123638"
2 Part 4
# Fit a GLS model with AR(1) correlation structure
gls_model <- gls(Sales ~ PDI, data = sales,
correlation = corAR1(form = ~ 1))
# Summary of the GLS model
summary(gls_model)
## Generalized least squares fit by REML
## Model: Sales ~ PDI
## Data: sales
## AIC BIC logLik
## 213.8913 220.4416 -102.9456
##
## Correlation Structure: AR(1)
## Formula: ~1
## Parameter estimate(s):
## Phi
## 0.05259568
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 12.426463 2.6742812 4.646655 0
## PDI 0.197734 0.0168734 11.718736 0
##
## Correlation:
## (Intr)
## PDI -0.982
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.7281062 -0.8876023 0.1800854 0.8982568 1.4696211
##
## Residual standard error: 3.027736
## Degrees of freedom: 40 total; 38 residual
# AIC and BIC comparison
cat("AIC for Cochrane-Orcutt Model:", AIC(sales_model1), "\n")
## AIC for Cochrane-Orcutt Model: 200.7121
cat("AIC for GLS Model:", AIC(gls_model), "\n")
## AIC for GLS Model: 213.8913
cat("BIC for Cochrane-Orcutt Model:", BIC(sales_model1), "\n")
## BIC for Cochrane-Orcutt Model: 205.7028
cat("BIC for GLS Model:", BIC(gls_model), "\n")
## BIC for GLS Model: 220.4416
"Since the Cochrane-Orcutt model has lower values for both criteria, it suggests that it
provides a better balance of goodness of fit and model simplicity for the data compared to
the GLS model.
"
## [1] "Since the Cochrane-Orcutt model has lower values for both criteria, it suggests that it\nprovides a better balance of goodness of fit and model simplicity for the data compared to\nthe GLS model.\n"
Q3 Part a)
load("/Users/fahadmehfooz/Downloads/tsa3.rda")
plot(jj, type="o", ylab="Quarterly Earnings per Share")
trend = time(jj) - 1970 # helps to 'center' time
Q = factor(cycle(jj) ) # make (Q)uarter factors
reg = lm(log(jj)~0 + trend + Q, na.action=NULL) # no intercept
model.matrix(reg)
## trend Q1 Q2 Q3 Q4
## 1 -10.00 1 0 0 0
## 2 -9.75 0 1 0 0
## 3 -9.50 0 0 1 0
## 4 -9.25 0 0 0 1
## 5 -9.00 1 0 0 0
## 6 -8.75 0 1 0 0
## 7 -8.50 0 0 1 0
## 8 -8.25 0 0 0 1
## 9 -8.00 1 0 0 0
## 10 -7.75 0 1 0 0
## 11 -7.50 0 0 1 0
## 12 -7.25 0 0 0 1
## 13 -7.00 1 0 0 0
## 14 -6.75 0 1 0 0
## 15 -6.50 0 0 1 0
## 16 -6.25 0 0 0 1
## 17 -6.00 1 0 0 0
## 18 -5.75 0 1 0 0
## 19 -5.50 0 0 1 0
## 20 -5.25 0 0 0 1
## 21 -5.00 1 0 0 0
## 22 -4.75 0 1 0 0
## 23 -4.50 0 0 1 0
## 24 -4.25 0 0 0 1
## 25 -4.00 1 0 0 0
## 26 -3.75 0 1 0 0
## 27 -3.50 0 0 1 0
## 28 -3.25 0 0 0 1
## 29 -3.00 1 0 0 0
## 30 -2.75 0 1 0 0
## 31 -2.50 0 0 1 0
## 32 -2.25 0 0 0 1
## 33 -2.00 1 0 0 0
## 34 -1.75 0 1 0 0
## 35 -1.50 0 0 1 0
## 36 -1.25 0 0 0 1
## 37 -1.00 1 0 0 0
## 38 -0.75 0 1 0 0
## 39 -0.50 0 0 1 0
## 40 -0.25 0 0 0 1
## 41 0.00 1 0 0 0
## 42 0.25 0 1 0 0
## 43 0.50 0 0 1 0
## 44 0.75 0 0 0 1
## 45 1.00 1 0 0 0
## 46 1.25 0 1 0 0
## 47 1.50 0 0 1 0
## 48 1.75 0 0 0 1
## 49 2.00 1 0 0 0
## 50 2.25 0 1 0 0
## 51 2.50 0 0 1 0
## 52 2.75 0 0 0 1
## 53 3.00 1 0 0 0
## 54 3.25 0 1 0 0
## 55 3.50 0 0 1 0
## 56 3.75 0 0 0 1
## 57 4.00 1 0 0 0
## 58 4.25 0 1 0 0
## 59 4.50 0 0 1 0
## 60 4.75 0 0 0 1
## 61 5.00 1 0 0 0
## 62 5.25 0 1 0 0
## 63 5.50 0 0 1 0
## 64 5.75 0 0 0 1
## 65 6.00 1 0 0 0
## 66 6.25 0 1 0 0
## 67 6.50 0 0 1 0
## 68 6.75 0 0 0 1
## 69 7.00 1 0 0 0
## 70 7.25 0 1 0 0
## 71 7.50 0 0 1 0
## 72 7.75 0 0 0 1
## 73 8.00 1 0 0 0
## 74 8.25 0 1 0 0
## 75 8.50 0 0 1 0
## 76 8.75 0 0 0 1
## 77 9.00 1 0 0 0
## 78 9.25 0 1 0 0
## 79 9.50 0 0 1 0
## 80 9.75 0 0 0 1
## 81 10.00 1 0 0 0
## 82 10.25 0 1 0 0
## 83 10.50 0 0 1 0
## 84 10.75 0 0 0 1
## attr(,"assign")
## [1] 1 2 2 2 2
## attr(,"contrasts")
## attr(,"contrasts")$Q
## [1] "contr.treatment"
Q3 Part b)
## log(jj)=β⋅trend+α1​Q1​(t)+α2​Q2​(t)+α3​Q3​(t)+α4​Q4​(t)+wt
"
β for the trend variable can be interpreted as the average annual
growth rate in the logged earnings per share. To interpret this in more understandable
terms
"
## [1] "\nβ for the trend variable can be interpreted as the average annual\ngrowth rate in the logged earnings per share. To interpret this in more understandable\nterms\n"
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
beta_trend <- coef(reg)["trend"]
beta_trend * 4
## trend
## 0.6686887
percentage_increase = beta_trend * 4 * 100
percentage_increase
## trend
## 66.86887
## average annual percentage increase in the earnings per share = 66.86887 %
Q3 Part c
coef_regs <- coef(reg)
coef_regs
## trend Q1 Q2 Q3 Q4
## 0.1671722 1.0527931 1.0809159 1.1510242 0.8822665
# percentage change from Q3 to Q4
percentage_change = ((coef_regs["Q4"] - coef_regs["Q3"])/coef_regs["Q3"]) * 100
percentage_change
## Q4
## -23.34944
print(paste("This corresponds to a percentage change of:", percentage_change, "%"))
## [1] "This corresponds to a percentage change of: -23.3494385958917 %"
"Reduced by 23.34%"
## [1] "Reduced by 23.34%"
Part d:
"
For model without Intercept:
Each alpha(coefficient) represents the mean log earnings for quarter i, directly, as there's no
baseline (intercept) to compare against.
The model assumes the mean log earnings can differ by quarter, independent of a base level.
For model with Intercept:
alpha0 (the intercept) represents the mean log earnings for a baseline quarter. If all
(t) indicators are 0 (which doesn't occur in this setup since one quarter must always be
present), it would theoretically represent the mean log earnings for this baseline.
However, in practice, since each time point is associated with a quarter,
alpha0 gets absorbed as part of the seasonal effects, making direct interpretation challenging.
Each alpha i now measures the difference in mean log earnings from the baseline quarter
represented by the intercept. This introduces redundancy because the seasonal effects are
fully captured by the quarters, making the intercept unnecessary or even misleading.
Problem with Including an Intercept:
Multicollinearity: Including an intercept with indicators for all quarters introduces
multicollinearity. The sum of the indicator variables for the quarters adds up to 1 for any
observation, making one of the indicators perfectly predictable from the others and the
intercept. This redundancy can lead to inflated standard errors and unstable estimates for the coefficients.
Interpretation Issues: With an intercept included, the coefficients for quarters are
interpreted relative to the baseline (the intercept), complicating the interpretation because
there is no non-quarter category — each time point falls into one of the four quarters.
"
## [1] "\nFor model without Intercept:\nEach alpha(coefficient) represents the mean log earnings for quarter i, directly, as there's no\nbaseline (intercept) to compare against.\nThe model assumes the mean log earnings can differ by quarter, independent of a base level.\n\nFor model with Intercept:\n\nalpha0 (the intercept) represents the mean log earnings for a baseline quarter. If all \n(t) indicators are 0 (which doesn't occur in this setup since one quarter must always be\npresent), it would theoretically represent the mean log earnings for this baseline.\nHowever, in practice, since each time point is associated with a quarter, \nalpha0 gets absorbed as part of the seasonal effects, making direct interpretation challenging.\n\n\nEach alpha i now measures the difference in mean log earnings from the baseline quarter\nrepresented by the intercept. This introduces redundancy because the seasonal effects are\nfully captured by the quarters, making the intercept unnecessary or even misleading.\n\nProblem with Including an Intercept:\n\nMulticollinearity: Including an intercept with indicators for all quarters introduces\nmulticollinearity. The sum of the indicator variables for the quarters adds up to 1 for any\nobservation, making one of the indicators perfectly predictable from the others and the\nintercept. This redundancy can lead to inflated standard errors and unstable estimates for the coefficients.\n\nInterpretation Issues: With an intercept included, the coefficients for quarters are \ninterpreted relative to the baseline (the intercept), complicating the interpretation because\nthere is no non-quarter category — each time point falls into one of the four quarters.\n"
Q3 Part e
# jj is our time series object and reg is the fitted model
plot(log(jj), type="o", ylab="Logged Earnings", main="Logged Earnings and Fitted Values")
lines(fitted(reg), col="red") # Adds fitted values in red
"This will plot the logged earnings per share and superimpose the fitted values from
your regression model.
"
## [1] "This will plot the logged earnings per share and superimpose the fitted values from\nyour regression model.\n"
plot(resid(reg), type="o", ylab="Residuals", main="Residuals of the Model")
abline(h=0, col="blue") # Adds a horizontal line at 0
"The variance in residual plot shows the errors are not constant, there is a lot of variance.
which means the model fit is not good. The error trend keeps going up and down. Still the residuals
are in a pretty small range for all the periods. Its not a perfect fit but still a good model."
## [1] "The variance in residual plot shows the errors are not constant, there is a lot of variance.\nwhich means the model fit is not good. The error trend keeps going up and down. Still the residuals\nare in a pretty small range for all the periods. Its not a perfect fit but still a good model."
acf(resid(reg), main="ACF of Residuals")
"ACF shows similar plot the model would have fitted well if the variance would
have been similar across the errors. It does not look like a white noise but still pretty close."
## [1] "ACF shows similar plot the model would have fitted well if the variance would\nhave been similar across the errors. It does not look like a white noise but still pretty close."
model_summary <- summary(reg)
# Extracting the adjusted R-squared value
adjusted_r_squared <- model_summary$adj.r.squared
print(adjusted_r_squared)
## [1] 0.9930649
"Adjusted r square is .99 which says its a good fit."
## [1] "Adjusted r square is .99 which says its a good fit."
Question 4
par(mfrow=c(3,1)) # plot the data
plot(cmort, main="Cardiovascular Mortality", xlab="", ylab="")
plot(tempr, main="Temperature", xlab="", ylab="")
plot(part, main="Particulates", xlab="", ylab="")
dev.new() # open a new graphic device
ts.plot(cmort,tempr,part, col=1:3) # all on same plot (not shown)
pairs(cbind(Mortality=cmort, Temperature=tempr, Particulates=part))
temp = tempr-mean(tempr) # center temperature
temp2 = temp^2
trend = time(cmort) # time
fit = lm(cmort~ trend + temp + temp2 + part, na.action=NULL)
summary(fit) # regression results
##
## 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(fit)) # ANOVA table (compare to next line)
## 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)))) # Table 2.1
## 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) # sample size
AIC(fit)/num - log(2*pi) # AIC
## [1] 4.721732
BIC(fit)/num - log(2*pi) # BIC
## [1] 4.771699
(AICc = log(sum(resid(fit)^2)/num) + (num+5)/(num-5-2)) # AICc
## [1] 4.722062
# Lag by 4 weeks prior
part_lagged <- c(rep(NA, 4), part[-((length(part)-3):length(part))])
fit_updated <- lm(cmort ~ trend + temp + temp2 + part + part_lagged)
summary(fit_updated)
##
## Call:
## lm(formula = cmort ~ trend + temp + temp2 + part + part_lagged)
##
## 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 ***
## part_lagged 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
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.608, Adjusted R-squared: 0.6041
## F-statistic: 154.5 on 5 and 498 DF, p-value: < 2.2e-16
Question 4 part B
# Calculate pairwise correlations
# Assuming cmort, trend are ts objects
# Converting ts objects to numeric vectors
cmort_vec <- as.numeric(cmort)
trend_vec <- as.numeric(trend)
# Combine into a data frame
data_df <- data.frame(
cmort = cmort_vec,
trend = trend_vec,
pt = part,
pt_lagged = part_lagged
)
pairs(data_df)
# Calculate and print pairwise correlations
cor_matrix <- cor(data_df, use="complete.obs") # Handles NA values if present
print(cor_matrix)
## cmort trend pt pt_lagged
## cmort 1.0000000 -0.45241813 0.44228958 0.52099930
## trend -0.4524181 1.00000000 -0.06784453 -0.08860669
## pt 0.4422896 -0.06784453 1.00000000 0.53405047
## pt_lagged 0.5209993 -0.08860669 0.53405047 1.00000000
"
Cmort and pt has a correlation of 0.44 and pt, pt-4 has a correlation of 0.53.
The relationship is stronger between pt and pt-4
"
## [1] "\nCmort and pt has a correlation of 0.44 and pt, pt-4 has a correlation of 0.53.\nThe relationship is stronger between pt and pt-4\n"