x = ts(scan("http://www.stat.ualberta.ca/~wiens/stat479/S&Sdatasets/clim-hyd.dat")) airtemp = ts(x[1:454]) dewpoint = ts(x[(454+1):(2*454)]) cloudcover = ts(x[(2*454+1):(3*454)]) windspeed = ts(x[(3*454+1):(4*454)]) precipitation = ts(x[(4*454+1):(5*454)]) inflow = ts(x[(5*454+1):(6*454)]) # Transform: inflow = log(inflow) precipitation = sqrt(precipitation) # Plot the data par(mfrow=c(3,2)) plot(airtemp) plot(dewpoint) plot(cloudcover) plot(windspeed) plot(precipitation) plot(inflow) title(main = "\n Shasta Lake inflow and predictors", outer = T) ###### Compute and plot marginal coherences: (compare with Figure 7.3) ###### win.graph() par(mfrow=c(3,2)) data = cbind(airtemp, dewpoint, cloudcover, windspeed, precipitation, inflow) datanames = c("airtemp", "dewpoint", "cloudcover", "windspeed", "precipitation", "inflow") L = 15 for (i in 1:5) { power = spec.pgram(data[,c(6,i)], kernel("daniell",(L-1)/2), plot = F) plot.spec.coherency(power, ci = 0, main = paste(datanames[6], 'vs', datanames[i])) f = qf(.999, 2, 2*(L-1)) c = f/(L-1+f) abline(h = c) } title(main = "\n Individual squared coherences with inflow", outer = T) ############################################################################## ############################################################################## source("http://www.stat.ualberta.ca/~wiens/stat679/R%20scripts/stoch.regr.R") ############################################################################## # This function computes the spectral matrix, F statistics and coherences, and plots them. # Returned as well are the coefficients in the impulse reponse function. # Enter, as the argument to this function, the full data matrix, and then the labels of the # columns of input series in the "full" and "reduced" regression models - enter NULL if there # are no inputs under the reduced model. # The response variable is the LAST column of the data matrix, and this need not be specified # among the inputs. Other inputs are alpha (test size), L (smoothing), M (number # of points in the discretization of the integral) and plot.which = "coh" # or "F", to plot either the coherences or the F statistics. ############################################################################## ############################################################################## ############################################################################## alpha = .01 L = 15 M = 30 data = cbind(airtemp, dewpoint, cloudcover, windspeed, precipitation, inflow) ####### F plots - compare with Figure 7.4 ######## win.graph() par(mfrow=c(2,2)) out.15 = stoch.regr(data, cols.full = c(1,5), cols.red = 5, alpha, L, M, plot.which = "F") title(main = "F to drop temp if precip remains", outer = F) out.25 = stoch.regr(data, cols.full = c(2,5), cols.red = 5, alpha, L, M, plot.which = "F") title(main = "F to drop dewpoint if precip remains", outer = F) out.35 = stoch.regr(data, cols.full = c(3,5), cols.red = 5, alpha, L, M, plot.which = "F") title(main = "F to drop cloudcover if precip remains", outer = F) out.45 = stoch.regr(data, cols.full = c(4,5), cols.red = 5, alpha, L, M, plot.which = "F") title(main = "F to drop windspeed if precip remains", outer = F) title(main = "\n F statistics", outer = T) ####### Coherency plots - compare the first with the top plot in Figure 7.5 ######## win.graph() par(mfrow=c(3,2)) coh.15 = stoch.regr(data, cols.full = c(1,5), cols.red = NULL, alpha, L, M, plot.which = "coh") title(main = "temp and precip vs. inflow", outer = F) coh.25 = stoch.regr(data, cols.full = c(2,5), cols.red = NULL, alpha, L, M, plot.which = "coh") title(main = "dewpoint and precip vs. inflow", outer = F) coh.35 = stoch.regr(data, cols.full = c(3,5), cols.red = NULL, alpha, L, M, plot.which = "coh") title(main = "cloudcover and precip vs. inflow", outer = F) coh.45 = stoch.regr(data, cols.full = c(4,5), cols.red = NULL, alpha, L, M, plot.which = "coh") title(main = "windspeed and precip vs. inflow", outer = F) coh.all = stoch.regr(data, cols.full = c(1:5), cols.red = NULL, alpha, L, M, plot.which = "coh") title(main = "all inputs vs. inflow", outer = F) title(main = "\n Multiple squared cohences", outer = T) ####### Impulse response functions - compare with the lower plots in Figure 7.5 ######## win.graph() par(mfrow=c(2,1)) S = seq(from = -M/2+1, to = M/2 - 1, length = M-1) plot(S, coh.15$Betahat[,1], type = "h", xlab = "", ylab = "temp", ylim = c(-.07, .07)) abline(h=0) plot(S, coh.15$Betahat[,2], type = "h", xlab = "", ylab = "precip", ylim = c(-.07, .07)) abline(h=0) title(main = "\n \n Impulse response functions relating inflow \n to temperature and precipitation", outer = T)