# The Frequency Domain Approach of a Time Series

t <- seq(0,7,length=100)
A <- 3
B <- -2
lam <- 0.2
y <- A*sin(2*pi*lam*t) + B*cos(2*pi*lam*t)
plot(t,y,type="l")

A1 <- 3
B1 <- -2
lam1 <- 0.2
A2 <- 1
B2 <- 3
lam2 <- 0.73
y <- A1*sin(2*pi*lam1*t) + B1*cos(2*pi*lam1*t) +
     A2*sin(2*pi*lam2*t) + B2*cos(2*pi*lam2*t)
plot(t,y,type="l")

# Example 3.1.1. (Star Data). To analyze physical properties of a pulsating
# star, the intensity of light emitted by this pulsar was recorded
# at midnight during 600 consecutive nights. The data are taken from
# Newton (1988).
star <- scan("star.txt")
star.160 <- star[1:160]
night <- 0:159
star.hat <- 17.07 - 1.86*cos(2*pi*(1/24)*night) + 6.82*sin(2*pi*(1/24)*night) + 
            6.09*cos(2*pi*(1/29)*night) + 8.01*sin(2*pi*(1/29)*night)
plot(night,star.160)
lines(night,star.hat)

star.data <- data.frame(star=star.160,cos24=cos(2*pi*(1/24)*night), sin24=sin(2*pi*(1/24)*night),
     cos29=cos(2*pi*(1/29)*night),sin29=sin(2*pi*(1/29)*night))
star.fit <- lm(star ~.,data=star.data)
plot(night,star.160)
lines(night,predict(star.fit))


# Priodogram and regression
n <- 10
T <- 2*pi*outer(1:n,(1:(n/2))/n)
X.cos <- cos(T)
X.sin <- sin(T)
round(t(X.cos)%*%X.cos,2)
round(t(X.sin)%*%X.sin,2)
round(t(X.cos)%*%X.sin,2)
X <- cbind(1,X.cos,X.sin[,-n/2])
round(t(X)%*%X,2)

# Plot 3.2.1a: Periodogram of the Star Data
n <- 600
T <- 2*pi*outer(1:n,(1:(n/2))/n)
X.cos <- cos(T)
X.sin <- sin(T)
X <- cbind(1,X.cos,X.sin[,-n/2])
b.hat <- as.vector(solve(t(X)%*%X)%*%t(X)%*%star)
plot((1:49)/600, b.hat[2:50]^2 + b.hat[2:50+n/2]^2,type="l")