\documentclass[12pt,titlepage]{article} \usepackage[dvips]{graphicx} \begin{document} \small %Spring 2006 %STAT 705 \begin{verbatim} B. Kedem, STAT 705, Fall Semester 2010, Time Series 2 a. Spectral Density Estimation for AR(1) Also: Cumulative Periodogram b. Random Telegraph Signal: Generation, ACF, Spectral Density Est. c. while() NOTE: Show example of raw periodogram. a. Spectral Density Estimation for AR(1) ========================================== A. Function to generate AR(1) ts of length n, parameter phi, and noise variance NoiseVar. The variance of the ts is NoiseVar/(1-phi^2). B. Autocorrelation Function of AR(1), phi=0.5. corr$acf[2] estimates phi. C. Spectral density of AR(1) with parameter phi, and noise variance NoiseVar. D. Estimated Spetral Density of AR(1). E. Cumulative Periodogram --------------------------------------- A. Function to generate AR(1) ts of length n, parameter phi, and noise variance NoiseVar. The variance of the ts is NoiseVar/(1-phi^2). ar1 <- function(n,NoiseVar,phi){ u <- sqrt(NoiseVar)*rnorm(n,0,1) x <- rep(NA,n) x[1] <- u[1] for(i in 2:n){x[i] <- phi*x[i-1] + u[i]} x} par(oma = c(1, 1, 5, 1)) par(mfrow = c(2,1), cex=0.8) z<- ar1(1000,1,0.5) ### AR(1): n=1000, NoiseVar=1, phi=0.5 ts.plot(z[1:300]) title("AR(1): n=300, NoiseVar=1, phi=0.5") z1<- ar1(1000,1,-0.5) ### AR(1): n=1000, NoiseVar=1, phi=-0.5 ts.plot(z1[1:300]) title("AR(1): n=300, NoiseVar=1, phi=-0.5") ##(FIG 1) mtext("ACF of AR(1), phi=0.5,-0.5", side = 3, outer = T, cex = 1.5, line=1) #OVERALL TITLE B. Autocorrelation Function of AR(1), phi=0.5 > corr <- acf(z) > names(corr) [1] "acf" "lag" "n.used" "type" "series" > corr$acf , , 1 [,1] [1,] 1.0000000000 <------ Rho(0) [2,] 0.5180638793 <------ Rho(1)=Least Squares Estimate of phi. [3,] 0.3067633050 [4,] 0.1933393272 [5,] 0.1627885289 [6,] 0.1432032082 [7,] 0.1096848343 [8,] 0.0639200109 [9,] 0.0078298110 [10,] -0.0066296156 [11,] -0.0320836670 [12,] -0.0446440149 [13,] -0.0229288801 [14,] -0.0096609308 [15,] -0.0053673497 [16,] -0.0114829070 [17,] -0.0356888854 [18,] -0.0450174614 [19,] -0.0564424616 [20,] -0.0388850345 [21,] -0.0035032387 [22,] 0.0326744560 [23,] 0.0382431709 [24,] -0.0008739067 [25,] -0.0258745207 [26,] -0.0121042661 [27,] -0.0182190777 [28,] 0.0324180664 [29,] 0.0176885699 [30,] 0.0189061701 [31,] -0.0141110941 > corr$acf[2] [1] 0.5180639 B'. Autocorrelation Function of AR(1), phi=-0.5 acf(z1) ##Of AR(1), phi=-0.5 C. Spectral density of AR(1) with parameter phi, and noise variance NoiseVar. f.ar1 <- function(NoiseVar,phi) ##AR(1) spectral density (NoiseVar/(2*pi))*(1-2*phi*cos(lambda)+phi^2)^(-1) EX 1: Spectral density of AR(1) with NoiseVar=1, phi=0.5 lambda <- seq(-pi,pi,length=400) ###For -pi < freq < pi plot(lambda,f.ar1(1,.5), type="l") EX 2: Spectral density of AR(1) with NoiseVar=1, phi=0.5 NoiseVar=3, phi=0.5 matplot(lambda,cbind(f.ar1(1,.5),f.ar1(3,.5)), type="l",lty=1:2) legend(-3.3,1.6,c("VarN=1","VarN=3"), lty=1:2,bty="n") #bty="n" cancels the frame mtext("Spectral Density: AR(1), phi=0.5", cex=1.2, line=1) EX 3: Spectral density of AR(1) with NoiseVar=1, phi=0.5 NoiseVar=1, phi=-0.5 matplot(lambda,cbind(f.ar1(1,.5),f.ar1(1,-0.5)), type="l",lty=1:2) legend(-2.3,0.6,c("phi=0.5","phi=-0.5"), lty=1:2,bty="n") #bty="n" cancels the frame mtext("Spectral Density: AR(1), phi=0.5,-0.5", cex=1.2, line=1) EX: 4: Plotting over the interval [0,pi] lambda <- seq(0,pi,length=400) ###For 0 < freq < pi matplot(lambda,cbind(f.ar1(1,.5),f.ar1(1,-0.5)), type="l",lty=1:2) legend(-2.3,0.6,c("phi=0.5","phi=-0.5"), lty=1:2,bty="n") #bty="n" cancels the frame mtext("Spectral Density: AR(1), phi=0.5,-0.5", cex=1.2, line=1) mtext("AR(1) Spectral densities in [-pi,pi, and [0,pi]", side = 3, outer = T, cex = 1.5, line=1) #OVERALL TITLE (FIG 2) D. Estimated Spetral Density of AR(1). Plotted on 10*log_10 dB scale in cycles per unit time over a freq. range [-0.5,0.5]. Thus the period is 1/freq. EX 1: Estimated Spetral Density using modified Daniell method with Sequence of moving average of 8,12,16 etc. Periodogram Ordinates. par(mfrow=c(2,2), oma = c(0, 0, 4, 0)) ##4 lines of outer margin at the top. frame() ##Clears the GUI ESPa <- spectrum(z,spans=c(3,5,7)) ##ESP=Estimated spectral density ESPb <- spectrum(z,spans=c(5,7,9)) ESPc <- spectrum(z,spans=c(9,11,17)) ESPd <- spectrum(z,spans=c(17,17)) mtext("Daniel: a. 3,5,7 b. 5,7,9 c. 9,11,17 d. 17,17", side = 3, outer = T, cex = 1.5, line=1) #OVERALL TITLE (FIG 3) > names(ESPa) [1] "freq" "spec" "coh" "phase" "spans" "filter" [7] "df" "bandwidth" "n.used" "series" "method" "taper" [13] "pad" "detrend" "demean" > length(ESPa$freq) [1] 500 > ESPa$freq[1:50] [1] 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.010 0.011 0.012 [13] 0.013 0.014 0.015 0.016 0.017 0.018 0.019 0.020 0.021 0.022 0.023 0.024 [25] 0.025 0.026 0.027 0.028 0.029 0.030 0.031 0.032 0.033 0.034 0.035 0.036 [37] 0.037 0.038 0.039 0.040 0.041 0.042 0.043 0.044 0.045 0.046 0.047 0.048 [49] 0.049 0.050 > ESPa$freq[451:500] [1] 0.451 0.452 0.453 0.454 0.455 0.456 0.457 0.458 0.459 0.460 0.461 0.462 [13] 0.463 0.464 0.465 0.466 0.467 0.468 0.469 0.470 0.471 0.472 0.473 0.474 [25] 0.475 0.476 0.477 0.478 0.479 0.480 0.481 0.482 0.483 0.484 0.485 0.486 [37] 0.487 0.488 0.489 0.490 0.491 0.492 0.493 0.494 0.495 0.496 0.497 0.498 [49] 0.499 0.500 > length(ESPa$spec) [1] 500 > ESPa$spec[1:50] [1] 3.720132 3.659553 3.597334 3.565221 3.542743 3.538213 3.632328 3.867884 [9] 4.146028 4.318514 4.282607 3.987530 3.460436 2.792231 2.118805 1.597785 [17] 1.356799 1.420830 1.660245 1.901774 2.051746 2.072970 1.953968 1.729719 [25] 1.500290 1.387820 1.461657 1.731698 2.173851 2.709244 3.200350 3.528170 [33] 3.659537 3.634342 3.514059 3.337990 3.131263 2.940865 2.801089 2.686121 [41] 2.555233 2.447937 2.530584 2.917244 3.489057 4.037049 4.426335 4.608851 [49] 4.568276 4.307746 > ESPa$spec[451:500] [1] 0.3980300 0.3611917 0.3412415 0.3375341 0.3557372 0.3979819 0.4490002 [8] 0.4920048 0.5222233 0.5419780 0.5550954 0.5648199 0.5695558 0.5629556 [15] 0.5481996 0.5375820 0.5330361 0.5299587 0.5346813 0.5594456 0.6047843 [22] 0.6582944 0.7079474 0.7429495 0.7529910 0.7340337 0.6867391 0.6217728 [29] 0.5563408 0.5015841 0.4636306 0.4429340 0.4385812 0.4465458 0.4518384 [36] 0.4461178 0.4378142 0.4380515 0.4495590 0.4689617 0.4948749 0.5250251 [43] 0.5506142 0.5595893 0.5441134 0.5053809 0.4523912 0.3973174 0.3529588 [50] 0.3351209 Spec. density on linear scale evaluated at frequencies in radians per unit time: par(mfrow=c(2,1)) z<- ar1(1000,1,0.5) ESPa <- spectrum(z,spans=c(3,5,7)) SpecDensity <- 10^(ESPa$spec/10) plot(2*pi*ESPa$freq, SpecDensity, type="l", main="Linear Scale in [0,pi]") mtext("Daniel(3,5,7): 10*Log_10 in [0,0.5] vs Linear Scale in [0,pi]", side = 3, outer = T, cex = 1.5, line=1) #OVERALL TITLE (FIG 4) z<- ar1(1000,1,0.5) ESPa <- spectrum(z,spans=c(3,3)) SpecDensity <- 10^(ESPa$spec/10) plot(2*pi*ESPa$freq, SpecDensity, type="l", main="Linear Scale in [0,pi]") mtext("Daniel(3,3): 10*Log_10 in [0,0.5] vs Linear Scale in [0,pi]", side = 3, outer = T, cex = 1.5, line=1) #OVERALL TITLE (FIG 4) E. Cumulative Periodogram ----------------------------- To plots a cumulative periodogram use: library(MASS) ###Must add. cpgram(ts, taper=0.1, main= paste("Series: ", deparse(substitute(ts)))) Useful for testing if the series is white noise: For WN the cumulative periodogram enters wn bounds. par(oma=c(0,0,4,0), mfrow=c(2,2)) z<- ar1(1000,1,0.8) ### AR(1): n=1000, NoiseVar=1, phi=0.8 cpgram(z,main="0.8") z<- ar1(1000,1,0.3) ### AR(1): n=1000, NoiseVar=1, phi=0.3 cpgram(z,main="0.3") z<- ar1(1000,1,0.0) ### AR(1): n=1000, NoiseVar=1, phi=0.0 cpgram(z,main="0.0") z<- ar1(1000,1,-0.8) ### AR(1): n=1000, NoiseVar=1, phi=-0.8 cpgram(z,main="-0.8") mtext("Cumulative Periodogram of AR(1)", side = 3, outer = T, cex = 1.5, line=1) #OVERALL TITLE (FIG 4) #Compare with ACF: par(oma=c(0,0,4,0), mfrow=c(2,2)) z<- ar1(1000,1,0.8) ### AR(1): n=1000, NoiseVar=1, phi=0.8 acf(z,main="0.8") z<- ar1(1000,1,0.3) ### AR(1): n=1000, NoiseVar=1, phi=0.3 acf(z,main="0.3") z<- ar1(1000,1,0.0) ### AR(1): n=1000, NoiseVar=1, phi=0.0 acf(z,main="0.0") z<- ar1(1000,1,-0.8) ### AR(1): n=1000, NoiseVar=1, phi=-0.8 acf(z,main="-0.8") mtext("Estimated ACF of AR(1)", side = 3, outer = T, cex = 1.5, line=1) #OVERALL TITLE (FIG 4) b. Random Telegraph Signal: Generation, ACF, Spectral Density Est. ================================================================== par(oma = c(0, 0, 4, 0)) ##4 lines of outer margin at the top par(mfrow=c(1,1)) ###Generate a Poisson process by generating the interarrival times which are independent Exp(lam). f <- function(a,x){a*exp(-a*x)} #Exponential(a) pdf #Histogram and pdf lam <- 0.1 x <- rexp(100,lam) ##Generated Exp(lam) hist(x,probability=T,cex=1.2,main="Exp(0.1)") q <- seq(0.001,60,length=500) lines(q,f(lam,q), type="l") mtext("Exponential Distribution", side = 3, outer = T, cex = 1.5, line=1) #OVERALL TITLE ###Sum: T1, T1+T2, T1+T2+T3,..... (same as cumsum(x)!!!) SUM <- 0 for(i in 1:length(x)){ SUM[i] <- sum(x[1:i])} ###Creating the indicator I[a < t <= b] NOTE: w <- function(x){(2 < x < 5)} does not work! w <- function(x){(2 < x ) & (x < 5)} ok! w1 <- function(a,b,t){ ##Indicator function 1*((a < t) & (t <= b))} #Sum of two indicators with different signs X<- 0 tpoints <- seq(0,SUM[50],length=5000) X <- w1(0,SUM[1],tpoints)-w1(SUM[1],SUM[2],tpoints) ###Sampling of Random Telegraph signal: Delta=1 > SUM[50] [1] 508.1134 <--- Get a ts of length > 500 tpoints <- seq(0,SUM[50],length=SUM[50]) X <- 0 #Summing all the indicators with alternating signs to get random telegraph square signal: for(i in 2:50){ X <- X + ((-1)^(i+1))*w1(SUM[i-1],SUM[i],tpoints)} X <- X+w1(0,SUM[1],tpoints) ts.plot(X,xlab="t", main="Random Telegraph: Lambda=0.1 Sampling Delta = 1") abline(0,0) acf(X) lines(tpoints,exp(-(2*lam*tpoints)),type="l") mtext("Random Telegraph: Lambda=0.1", side = 3, outer = T, cex = 1.5, line=1) mtext("Estimated vs Theoretical ACF", side = 3, outer = T, cex = 1.5, line=-4) #Spec. density of Rndom Telegraph on linear scale evaluated at frequencies in radians per unit time. Lambda=0.1, Dela t = 1. par(oma = c(0, 0, 4, 0)) par(mfrow=c(2,1)) ESPa <- spectrum(X,spans=c(3,5,7)) #First in Hz on [0,0.5] #Now in radians/unit time SpecDensity <- 10^(ESPa$spec/10) plot(2*pi*ESPa$freq, SpecDensity, type="l", main="Linear Scale in [0,pi]", cex=1.3,ylab="") mtext("Random Telegraph: Lambda=0.1, Delta t =1 Daniel(3,5,7): 10*Log_10 in [0,0.5] vs Linear Scale in [0,pi]", side = 3, outer = T, cex = 1.5, line=1) #OVERALL TITLE ##More refined plot: Sample often. tpoints <- seq(0,SUM[50],length=5000) X <- 0 #Summing all the indicators with alternating signs to get random telegraph signal: dev.off() for(i in 2:50){ X <- X + ((-1)^(i+1))*w1(SUM[i-1],SUM[i],tpoints)} X <- X+w1(0,SUM[1],tpoints) ts.plot(X,xlab="t", main="Random Telegraph: Lambda=0.1 Sampling Delta = 0.085") abline(0,0) #Better plot: time sxis referes to actual time plot(tpoints,X,type="l",xlab="t",cex=1.3,xlim=c(0,450)) abline(0) mtext("Random Telegraph: Lambda=0.1, Delta t = 0.085 ", side = 3, outer = T, cex = 1.5, line=-4) acf(X) ##Simpler way to get plot of random telegraph signal using type="s" or "S" for stairsteps ---------- lam <- 0.1 x <- rexp(50,lam) SUM <- cumsum(x) plot(SUM,rep(c(-1,1),25),type="s",xlab="TIME",ylab="") #Lower case "s" plot(SUM,rep(c(-1,1),25),type="S",xlab="TIME",ylab="") #Cap "S" Compare with: plot(cumsum(x),rnorm(50),type="s",xlab="TIME",ylab="") c. while() =========================== n <- 5 while(n > 0) { cat("n is still greater than 0\n") n <- n - 1 } n is still greater than 0 n is still greater than 0 n is still greater than 0 n is still greater than 0 n is still greater than 0 \end{verbatim} \end{document}