Analysis of Spatio Temporal Data : Lab2

Avipsa Roy, MSc Geoinformatics : WISE 2015/16

1. Compute the periodogram of the raw co2 time series.

N = length(co2)
I = abs(fft(co2)/sqrt(N))^2
## Scale periodogram
P = (4/N) * I  ## scaled periodogram
f = (0:floor(N/2))/N
plot(f, I[1:((N/2) + 1)], type = "o", xlab = "frequency", ylab = "", main = "Raw periodogram without trend removal")

spec.pgram(co2, detrend = FALSE, ci=0, col = "blue", main = "without (blue) and with (red) trend removal")
spec.pgram(co2, detrend = TRUE, ci = 0, col = "red", add = TRUE)

2. Use function stl() to remove the trend/seasonal component from the CO2 data. Recompute the periodogram and compare it to that in the previous exercise.

model = stl(co2, s.window = "periodic")
trend = model$time.series[, "trend"]
co2.no.trend = co2 - trend
plot(co2.no.trend)

N = length(co2.no.trend)
I = abs(fft(co2.no.trend)/sqrt(N))^2
P = (4/N) * I
f = (0:floor(N/2))/N
plot(f, I[1:((N/2) + 1)], type = "h", xlab = "frequency", ylab = "", main = "Periodogram of CO2 data after trend removal", col = "blue")

3. Draw 500 samples from the standard normal distribution. Compute the peridogram of this white noise time series.

noise = rnorm(1:500)
N = 500
I = abs(fft(noise)/sqrt(N))^2
P = (4/N) * I
f = (0:floor(N/2))/N
plot(f, I[1:((N/2) + 1)], type = "o", xlab = "frequency", ylab = "", main = "Periodogram of White Noise")

4. Simulate 500 time steps from the following cosine model: Yt=A⋅cos(2πωt+ϕ)+et with ω=1/50, A=2, ϕ=0.6π and et∼N(0,σ2) for σ2=1 and σ2=5 each. Compute the periodograms and compare them also with the above periodogram.

noise = rnorm(500, 0, sqrt(5))
t = 1:500
et1 = seq(0,1,0.002)
et2 = seq(0,5,0.01)
model1 = 2 * cos(1/50 * 2 * pi * t + 0.6 * pi) + et1[1:500]
model2 = 2 * cos(1/50 * 2 * pi * t + 0.6 * pi) + et2[1:500]
ysim1 = model1 + noise
ysim2 = model2 + noise
plot(ysim1, main = "Periodogram for cosine model with sigma^2=1")
lines(model1, type = "l", col = "red", lwd = 2)  

plot(ysim2, main = "Periodogram for cosine model with sigma^2=5")
lines(model2, type = "l", col = "blue", lwd = 2)

5. Simulate from an AR(1)-model using α=0.8 and α=−0.8. Try to explain the differences in the periodograms of these two series (R function arima.sim()). Try function spec.gram() to compute the periodogram.

ar1 = arima.sim(n = 500, list(ar = 0.8))
plot(ar2, type = "l", col = "blue", main = "AR(1) model with alpha=0.8")

ar2 = arima.sim(n = 500, list(ar = -0.8))
plot(ar2, type = "l", col = "red", main = "AR(1) model with alpha=-0.8")

spec.pgram(ar1, col = "blue", ci = 0, type = "s", main = "Raw periodogram of AR(1) model with alpha=0.8 and alpha = -0.8")
spec.pgram(ar2, col = "red", ci = 0, type = "s", add = TRUE)

6. Read more about leakage and antialiasing. Write down a brief explanation of these phenomena.

7. Use function spectrum as an alternative to estimate the contributing frequencies of the above time series.