#HANDS-ON INTERMEDIATE ECONOMETRICS USING R #Second Edition, 2022, World Scientific Publishers #by Hrishikesh D Vinod #Following snippets are from Chapter 2 #If you use them, please give credit to Prof. Vinod #R2.1.1 Ejemplos: GDP, precio IBM, retornos precios S&P 500 install.packages("Rtools") install.packages("quantmod") library(quantmod) install.packages("getSymbols") install.packages("tseries") library(tseries) quantmod::getSymbols("GDPC1",src = "FRED") # start=(c(1947,1)) una columna de serie rgdp=ts(GDPC1,start=1947,freq=4) # frecuencia trimestral en GDPC1 ibm <- get.hist.quote("ibm", quote="Adj", start="1970-01-01", retclass="zoo", compression="m") # permite bajar datos de Internet sp500 <- get.hist.quote("^GSPC", quote="Adj", start="1970-01-01", retclass="zoo", compression="m") # permite bajar datos de Internet ribm <- 100*diff(log(ibm)) # retornos par(mfrow=c(2,2)) # crea un marco gráfico de 2 por 2 (cuatro gráficos) plot(rgdp, main="GDP real en US", ylab="GDP real") plot(ibm, main= "Precios mensuales de las acciones de IBM") plot(ribm, main="Retornos mesuales de las acciones de IBM") plot(sp500, main="Precios mensuales de las acciones de S&P 500") #R2.1.1b Implementación de Box Cox. Production #function example of boxcox() is in Chapter 11 snippets #R11.2.1############## #and #R11.2.2. set.seed(22) y=sample(1:10) x1=sample(1:10) x2=sample(11:20) reg=lm(y~x1+x2) install.packages("EnvStats") library(EnvStats) EnvStats::boxcox(reg, objective.name = "Shapiro-Wilk",optimize = TRUE) MASS::boxcox(Volume ~ log(Height) + log(Girth), data = trees,lambda = seq(-0.25, 0.25, length = 10)) #R2.1.2 Descomposición estacional de los datos de pasajeros de aerolíneas summary(AirPassengers) # datos iniciales dair <- decompose(AirPassengers) # objeto R tiene descomposición plot(dair) # produce 4 gráficos de los datos originales sum <- dair$seasonal+dair$trend+dair$random cbind(sum,AirPassengers, dair$seasonal, dair$trend, dair$random) #crea matriz sair <- stl(AirPassengers,s.window="periodic") #descompone una serie en estacional #tendencia y componentes irregulares plot(sair) # gráfica más refinada, evita valores perdidos y ajustes polinomiales mtx <- as.matrix(sair$time.series) dimnames(mtx) sum2 <- apply(mtx,1,sum) # computa suma de mtx en una columna airseas <- mtx[,"seasonal"] # nuevo vector columna para futuro uso en ajuste estacional airtrend<- mtx[,"trend"] # nuevo vector columna para futuro uso en ajuste tendencia cbind(airseas,airtrend,mtx[,"remainder"],sum2,AirPassengers ) # une columnas ########################################################################### #R2.1.3 Filtro de Hodrick-Prescott data de Pasajeros AirPassenger NO HAY PERMISO DE ACCESO #require(graphics) #devtools::install_github("chenyang45/BoostedHP/BoostedHP") #bdair = BoostedHP(AirPassengers, lambda = 1600*3^4) #par(mfrow=c(3,1)) #plot(bdair$raw_data, ylab="data") #plot(bdair$cycle, ylab="cycle") #plot(bdair$trend,ylab="trend") ############################################################################ #R2.3.1 Gráfico seno coseno, dos ciclos t=1:100; thet=2*pi*(t-1)/100 plot(cos(thet),typ="l", main="Ondas básicas seno coseno") lines(sin(thet),typ="l", lty=2) thet=2*pi*(t-1)/50 #halving denominator doubles Cycle No. plot(cos(thet),typ="l", main="Dos ciclos de onda seno coseno") lines(sin(thet),typ="l", lty=2) #R2.3.2 Función en R para crear un AR(2) con raíces complejas desde N(0,1) ar2.cmplex=function(phi1, phi2,n=100){ y=rep(NA,n);#place to store series set.seed(239); rn=rnorm(n) y[1]=rn[1] y[2]=phi1*y[1]+rn[2] for (t in 3:n) { y[t]=phi1*y[t-1]+phi2*y[t-2]+rn[3] } #loop ends here return(y)} #end of the function plot(ts(ar2.cmplex(0.2,-.97)),main="Sinusoide amortiguada de AR(2) con raíces complejas") # esto entrega un ciclo amortiguado. Note que si no tenemos # phi2<0, no podemos tener ciclos. Además, necesitamos # phi1^2 < 4*|phi2| para que el radical sea negativo #R2.4.1 Función de autocorrelación (ACF) para el GNP real library(tseries) data(NelPlo)# Nelson-Plosser Macroeconomic Time Series summary(NelPlo) acf(gnp.real) #R2.4.2 Una función en R modifica la acf la función disponible acf2=function(data, maxlag=NULL){ num=length(data) if (is.null(maxlag)) maxlag=ceiling(10+sqrt(num)) mp1=maxlag+1 ACF=acf(data, maxlag, plot=F)$acf[2:mp1] PACF=pacf(data, maxlag, plot=F)$acf LAG=1:maxlag minA=min(ACF) minP=min(PACF) U=2/sqrt(num) L=-U minu=min(minA,minP,L)-.01 par(mfrow=c(2,1), mar = c(3,3,2,0.8), oma = c(1,1.2,1,1), mgp = c(1.5,0.6,0)) plot(LAG, ACF, type="h",ylim=c(minu,1), main=paste("ACF, PACF for ",deparse(substitute(data)))) abline(h=c(0,L,U), lty=c(1,2,2), col=c(1,4,4)) plot(LAG, PACF, type="h",ylim=c(minu,1)) abline(h=c(0,L,U), lty=c(1,2,2), col=c(1,4,4)) return(cbind(LAG, ACF, PACF))}#función termina aquí #Corra #R2.1.2 y renombre airseas seasonal.AirPassengers <- airseas acf2(seasonal.AirPassengers) # PACF para rezagos 1 y 2 están # fuera de bandas sugiere AR(2). ACF todos los rezagos están afuera. #R2.4.3 Ejemplos de simulaciones AR y MA. Necesita correr #R2.4.2 AR0.7 <- arima.sim(list(order = c(1,0,0), ar = 0.7), n = 50) summary(AR0.7) #c(1,0,0) 1 es AR(1), 0 es no diferencias y 0 es no media móvil ts.plot(AR0.7) acf(AR0.7) #función de autocorrelación declina suavemente pacf(AR0.7) # para pacf ar(1) es el model sugerido # copie codigo acf2 de #R2.4.2 acf2(AR0.7)# sólo un rezago es significativo AR(1) # Ejemplo de una simulación de media móvil set.seed(345) MA0.7 <- arima.sim(list(order = c(0,0,1), ma = 0.7), n = 50) ts.plot(MA0.7) # debe tener R2.4.2 en memoria acf2(MA0.7)#sugiere MA(1) ya que sólo un rezago de ACF es significativo. #R2.4.4 Componentes estacionales para la data de Airlines #debe tener R2.4.2 en memoria acf2(airseas) #sugiere AR(2) PACF tiene 2 rezagos acf2(airtrend) #suggests AR(1), PACF tiene 1 rezago #R2.4.5 Adecuando un modelo ARIMA para el GNP real desestacionalizado library(tseries) data(NelPlo) # Nelson\(-\)Plosser Macroeconomic Time Series summary(NelPlo) plot(gnp.real, main="US Real GNP") arima(log(gnp.real), order=c(2,0,0)) #OUTPT: Coefficients and standard errors for Real GNP: # ar1 ar2 intercept # 1.4377 -0.4408 1.7595 #s.e. 0.1001 0.1006 0.1665 # Fitting ARIMA to deseasonal parts of airline passenger data arima(airseas, order=c(2,0,0)) #arma needs library terseries #coefficients phi1 and phi2 suggest imaginary roots in #seasonal part arima(airtrend,c(1,0,0)) #coefficient phi exceeds 1, #nonstationary, #need differencing with d=1 to make it stationary arima(airtrend, c(1,1,0)) #set d=1 to do differencing #R2.4.6 Ejemplos simulados de ARIMA con diagnóstico MA0.7 <- arima.sim(list(order = c(0,0,1), ma = 0.7), n = 50) fit2 = arima(MA0.7, order=c(0,0,1));fit2 tsdiag(fit2) # dibujamos 3 cuadros #model fit2 #OUTPUT: Estimated MA(1) coef is close to 0.7 as expected # ma1 intercept # 0.6532 -0.3369 #s.e. 0.1216 0.2215 #Intercept is not significant #sigma^2 estimate= 0.9115,log likelihood = -68.91,aic = 143.82 #R2.4.6b ARMA estimation and diagnostics: AR simulation set.seed(345) AR0.7 <- arima.sim(list(order = c(1,0,0), ar = 0.7), n = 50) fit = arima(AR0.7, order=c(1,0,0)) tsdiag(fit)#3 cuadros muestran que el modelo correcto #de residuos estandarizados no muestra un particular patrón #ACF de residuos muestran insignificantes autocorrelaciones de errores #p-values for Ljung-Box estadísticos son mayores que 0.05. De esta forma el #producto es similar a un MA(1). Resultado omitido por brevedad. #R2.4.7 Ejemplo 2 ARIMA(p,d,q)x(P,D,Q) para datos AirPassengers a1 <- arima(AirPassengers, order = c(1, 1, 0), #order d=1=D seasonal=list(order = c(2, 1, 0), period = 12)); a1 tsdiag(a1) #modelo a1: d=1=D,p=1,P=2 es bueno a2 <- arima(AirPassengers, order = c(1, 0, 0), #order d=0,D=1 seasonal=list(order = c(2, 1, 0), period = 12)); a2 tsdiag(a2)#modelo a2 con d=0, D=1, P=2 es rechazado a3 <- arima(AirPassengers, order = c(1, 1, 0), seasonal=list(order = c(2, 0, 0), period = 12)); a3 tsdiag(a3) #a3 falla, no estacionariedad estacional #R2.4.8 Algunos test: AIC, Jarque-Bera, Ljung-Box #Los resultados no recomiendan modelo “a2” incluso si tiene #mejores propiedades AIC, BIC y normalidad que “a1”. La prueba gráfica de #normalidad de los residuos de regresión, a gráfica cuantil-cuantil #frente a la normal estándar mediante el comando “qqnorma”. #Por razones de brevedad, omitimos muchos de esos detalles. #Asuma que R2.4.7 incluyendo a1 y a2 están en memoria #tsdiag(a2) AIC(a1);AIC(a2) #imprime valores del criterio de información de Akaike #minimizar criterio AIC es recomendado. AIC(a1, k = log(length(AirPassengers))) AIC(a2, k = log(length(AirPassengers))) #sugiere modelo a2 #Diagnóstico: Jarque\(-\)Bera test de normalidad library(tseries) jarque.bera.test(residuals(a1))# p-value <0.05, #no son normal jarque.bera.test(residuals(a2))# p-value <0.05, #no son normal #Ljung-Box tests with p-values. for (ai in 1:2){ if (ai==1) aa=a1 if (ai==2) aa=a2 Mjb <- matrix(nrow=24, ncol=2) for(i in 1:24){ jb <- Box.test(residuals(aa), lag=i, type="Ljung-Box") Mjb[i,1] <- jb$statistic Mjb[i,2] <- jb$p.value } Mjb <- data.frame(Mjb) dimnames(Mjb) <- list(paste("Lag", 1:24, sep="."), c("stat.", "pval")) print(c("model= a",ai), quote=F) print(head(Mjb)) } #claramente sugiere modelo a1 como normal según Ljung-Box #R2.5.1 Gráficos simulados de ruido blanco set.seed(3467)#semilla rn=rnorm(100)#Crea la serie de ruido blanco con 100 observaciones plot.ts(rn, main="Serie de tiempo de ruido blanco") #R2.5.2 función de autocorrelación de data airline acf(AirPassengers, type="covariance", plot=F) #plot=F necesario para ver autocovarianzas #reporta la primera dos covarianza como 14291 y 13549 #R2.5.3 Camino aleatorio simulado set.seed(349); rn=rnorm(100); rwalk=cumsum(rn)#cumulative sum plot(rwalk) #R2.5.4 Gráfico de camino aleatorio simulado con tendencia Fig. 2.12 set.seed(349); rn=rnorm(100); rwalk=cumsum(rn)#suma acumulada set.seed(349); rn=rnorm(100); mu=rep(0.5,100) drift.rwalk=cumsum((rn+mu))#suma acumulada par(mfrow=c(2,1)) ; plot(rwalk, typ="l", main="Camino aleatorio") plot(drift.rwalk, typ="l", main="Camino aletorio con Drift =0.5") #R2.5.5 Tendencia polinomial t=seq(1,100,1)# secuencia de números 1 a 100 set.seed(3467)#permite replicación rn=rnorm(100)# creates a vector of 100 unit normals #Next define 3 vectors of order p polynomial data in time t: y1=1+2*t+rn; #(equation (2) with p=1 here) y2=1+2*t+3*t^2+rn; # (p=2 regressors here) y3=1+2*t+3*t^2+4*t^3+rn;# (p=3 regressors here) #Compute the first difference of y1 plot(y1, typ="l")#draw line plot of the linear trend fdy1=diff(y1,1) #construct a series for first difference of y1 lines(fdy1) # note first differenced series removes trend # removing trend means flattening the trend line. fdy2=diff(y2,2) fdy3=diff(y3,3)#series for third difference of cubic #polynomial data y3 # Plot all together with headings par(mfrow=c(3,1)) plot(y1, typ="l", main="Linear Trend and first differences(dashed line)") lines(fdy1, lty=2) # plot(y2, typ="l", main="Quadratic Trend and second differences(dashed line)") lines(fdy2, lty=2) # plot(y3, typ="l", main="Cubic Trend and third differences (dashed line)") lines(fdy3, lty=2) # #R2.5.6 Explosive growth trend t=1:100; yt=1.001^(0.1+0.5*t+0.5*t^2) plot(yt, typ="l", main="explosively growing trend") #R2.5.7 S-shaped trend ft = function(b0,b1,t) {ft = 1 / (1+exp( b0 - b1*t))} t=seq(-5,5,.1); ftt=ft(t,b0=0.5,b1=1); plot(ftt,t, main="Three S-shaped trends", typ="l") #solid line lines(ft(t,b0=0.5,b1=0.5),t, lty=2)#dashed line lines(ft(t,b0=0.5,b1=-0.5),t, col="blue", lty=3)#dotted line #R2.5.8 Understanding Amplitude and Phase angle #Imagine 5 years of monthly data (freq=12) yrs=5;freq=12; beta=0.5; psi=pi/3 #choose parameters t=1:(freq*yrs) tm=(t-1)/freq thet=2*pi*tm thet2=2*pi*tm + psi plot(cos(thet),typ="l",main="Original vs halved amplitude and changed phase") lines(beta*cos(thet2),typ="l", lty=2) #plots second wave #R2.5.9 Fitting Sinusoidal Trends to Data air=AirPassengers #use some data always available in R yrs=12;freq=12 t=1:(freq*yrs) tm=(t-1)/freq tpft=2*pi*tm coss=cos(tpft) #create first regressor in (2.5.13) sinn=sin(tpft)#create second regressor in (2.5.13) reg1=lm(air~coss+sinn) #Create regr. object reg1 fitt=fitted(reg1)#fitted values from regression matplot(cbind(air,fitt), typ="l") #plot data and fitted values # see that this regression does not allow for upward trend r1=resid(reg1) plot(resid(reg1),typ="l") # residuals have not detrended well # try another regression reg2=lm(r1~t) #new regression output object called reg2 plot(resid(reg2), typ="l")#Detrending removed up trend #following uses seasonal decomposition before trend fitting dct=decompose(air)$trend reg3=lm(dct~t)#new regression output object called reg3 lines(resid(reg3),typ="l");AIC(reg2);AIC(reg3) #Now compare reg2 with reg3 (Reject reg3) using AIC AIC(reg2); AIC(reg3) #minimize AIC by reg3. Box.test(resid(reg2), type="Ljung")#reject indep. errors Box.test(resid(reg3), type="Ljung")#reject indep. errors #R2.5.10 Plotting trend stationary and difference stationary #data t=1:100#a sequence of numbers from 1 to 100 set.seed(256)#fixed seed helps replication rn=rnorm(100)# create a vector of 100 unit normals #Next define vectors of linear trend data y1=ts(2*t+rn) #driftless trend as time series ts reg1=lm(y1~t) r1=ts(resid(reg1)) # detrending by subtraction of fitted trend= #residuals. ts makes them time series, easy to intersect dy=diff(y1,1);tsd=ts.intersect(r1,dy) matplot(tsd, typ="l", main="LTSP (residual) & FDSP (differenced as dashed line)",ylab="resid & differenced") y1b=ts(1+t+rn)# linear trend with drift reg1=lm(y1b~t-1) r1=ts(resid(reg1)) #drifted linear trend removed dy=diff(y1b,2);tsd=ts.intersect(r1,dy) matplot(tsd, typ="l", main="LTSP (drifted trend) & FDSP (differenced=dashed)",ylab="resid & differenced") # solid curve (LTSP) dashed (FDSP) both seem stationary #R2.5.11 Unit root testing X=AirPassengers # use your own data here library(urca);library(tseries);library(fUnitRoots) ur.df(X) # this gives the test statistics but not P values #this problem can be solved by calling additional package adf.test(X) adfTest(X) # accepts the null of unit root gives P value urdfTest(X) #gives greater details and graphs for ADF tests urersTest(X)# Elliott-Rothenberg-Stock tests unit root null urkpssTest(X) #KPSS test with stationarity as the null # Phillips and Perron test use following urppTest(X)#does beautiful graphs in addition to testing urzaTest(X) # For Zivot and Andrews test with figures urspTest(X) # For Schmidt-Phillips test with figures #R2.6.1 Mean Reversion in GNP library(urca);data(npext);attach(npext);summary(gnpperca); n=length(gnpperca);gnp=na.omit(gnpperca) bigt=length(gnp); bigk=bigt/2; vhat=rep(NA,bigk); se=rep(NA,bigk) for (k in 1:bigk) { varnum=var(diff(gnp, k)) varden= var(diff(gnp, 1)) # now implement Eq. (2.6.1) vhat[k] =(varnum*bigt)/ (k*varden*(bigt-k+1) ) se[k]=vhat[k]*sqrt( (4*k) / (3*bigt))} vhat; se vlow=vhat-se #compute lower conf limit vup=vhat+se #compute upper conf limit m1=min(vlow)#find limits for plotting m2=max(vup) plot.ts(vhat,ylim=c(m1,m2), main="Nonparametric Mean Reversion in per capita GNP Variance Ratio decreases as k increases", xlab="Lag, k", ylab="Variance ratio and confidence limits") lines(vlow, lty=2)#draw line lower limit lines(vup, lty=3) ac=acf(gnp, plot=F)#acf details without plots ac; names(ac); min(ac$acf) #extract min rho-k #R2.7.1 Plotting flat power spectrum set.seed(2345); x=rnorm(100) spec.ar(x)# this gives a flat power spectrum assuming AR model spectrum(x, spans=3) # gives sample periodogram which is seen #to lie inside the conf interval, flat spectrum near zero. #R2.7.2 Plotting periodogram and spectrum set.seed(345) AR0.7 <- arima.sim(list(order = c(1,0,0), ar = 0.7), n = 50) summary(AR0.7) par(mfrow=c(3,1)) #range (-frequency(x)/2, +frequency(x)/2] spectrum(AR0.7) spectrum(AR0.7, spans=9) spec.ar(AR0.7, main="Power spectrum: AR(1) with coefficient= 0.7") #R2.8.1 Fitting GARCH to European stock data and IBM library(tseries);data(EuStockMarkets) dax <- diff(log(EuStockMarkets))[,"DAX"] dax.garch <- garch(dax) # Fit a GARCH(1,1) to DAX returns summary(dax.garch) # ARCH effects are filtered. However, plot(dax.garch) # conditional normality seems violated ibm <- get.hist.quote("ibm", quote="Adj", start="1970-01-01", retclass="zoo", compression="m") plot(ibm, main="monthly prices of IBM stock") # Now we can compute monthly returns ribm <- 100*diff(log(ibm)) Ribm=as.numeric(ribm)#without this, next line gives NAs ribm.garch <- garch(Ribm, order = c(1,1))# fit GARCH(1,1) summary(ribm.garch) plot(ribm.garch) # conditional normality seems violated qqnorm(resid(ribm.garch)) r2=residuals(ribm.garch) r3=r2[r2==NA]; r4=setdiff(r2,r3) #remove missing residuals jarque.bera.test(r4)# p-value near zero, Reject normal resid library(fGarch) #load this package in R session. fit=garchFit(~garch(1,1),data=ribm) plot(fit) #gives 10 types of sophisticated plot options library(rugarch); data(dmbp); rets <- ts(dmbp$V1) spec=ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1), variance.targeting = FALSE), mean.model = list(armaOrder = c(1, 1), include.mean = TRUE, archm = FALSE, archpow = 1, arfima = FALSE), distribution.model = "norm") fit=ugarchfit(spec=spec, data=rets);fit plot(fit,which="all")#plots 12 diagnostics #R2.9.1 BDS & Teraesvirta Neural Network test library(fNonlinear); x=AirPassengers bdsTest(x); #Since p-values are almost zero, reject iid assump. for x data set.seed(4711); x = tentSim(n = 100) plot(x, main = "Tent Map Simulation") bdsTest(x, m = 3); tnnTest(x) set.seed(459)#create mixture density in R x = c(rnorm(50), runif(50)) plot(x, type = "l", main = "Mixture of Normal and Uniform Non-iid Time Series") bdsTest(x, m = 3);tnnTest(x) #R2.9.3 Structural Shift UK examples. library(strucchange) data("UKDriverDeaths") seatbelt <- log10(UKDriverDeaths) seatbelt <- cbind(seatbelt, lag(seatbelt, k = -1), lag(seatbelt, k = -12)) colnames(seatbelt) <- c("y", "ylag1", "ylag12") seatbelt <- window(seatbelt, start = c(1970, 1), end = c(1984,12)) plot(seatbelt[,"y"], ylab = expression(log[10] (casualties)), main="Structural shift in UK auto casualties") #R2.10.1 Nile river flows have long memory arima(Nile, order=c(1,0,0)) #yields AR(1) coefficient of #0.5063 kk=1:20; NileAR1=0.5063^kk; plot.ts(NileAR1) par(mfrow=c(2,1)) acf(Nile, main="Autocorrelations for Nile river flows") plot.ts(NileAR1, main="Autocorrelations for AR(1) model #for Nile river flows") #R2.10.2 Fractional diff estimate Nile flows and GNP library(fracdiff) fdGPH(Nile, bandw.exp = 0.5)# GPH estimate d=0.39 fr1=fracdiff(Nile) summary(fr1)#gives the ML estimate d=0.3639127 confint(fr1) #95% conf int.=(0.3638995, 0.3639258) library(tseries); data(USeconomic) summary(USeconomic) dgnp=diff(log(GNP),1) fracdiff(log(GNP))# d almost 0.5 is unstable level dgnp=diff(log(GNP),1)#compute first difference fracdiff(dgnp)# this d =0.2233329 is more plausible #R2.11.1 Predicting ARIMA(p,d,q)x(P,D,Q)for seasonal a1=arima(AirPassengers, order = c(1, 1, 0), seasonal = list(order = c(2, 1, 0), period = 12)) predict(a1, n.ahead=6) #R2.11.2 Example of forecasting US accidental deaths library(forecast) par(mfrow = c(4, 1)) #set up 4 plots in one plot(USAccDeaths, main="US accidental deaths") #plot original series accd.fcast <- meanf(USAccDeaths, h=10, fan=T) #h=number of periods plot(accd.fcast) # random walk forecast rwf NEXT plot(rwf(USAccDeaths, h=10, drift=FALSE, level=c(80,95), fan=T)) fit <- ets(USAccDeaths) #exponential smoothing state-space model above plot(forecast(fit)) #R2.11.3 VAR forecasting Canadian macroeconomic data library(vars); data(Canada)#load macroeconomic data for Canada par(mfrow = c(4, 1)) #set up 4 plots in one #prod= labor productivity; e=employment; U=unemployment rate # rw = real wage VAR vector autoregression estimated as var.2c <- VAR(Canada, p = 2, type = "const") var.2c.prd <- predict(var.2c, n.ahead = 8, ci = 0.95) fanchart(var.2c.prd) # does beautiful fanplot for forecasting #R2.12.1 AR(2) polynomial with complex roots y=c(10, 8, 22, 17, 21, 28, 10, 23, 24, 9) a2=arima(y, order = c(2, 0, 0)) #AR order p=2, differencing order d=0, and MA order q=0 summary(a2) a2$coef[1:2] phi=as.numeric(a2$coef[1:2]) rad=phi[1]^2 + 4*phi[2] rad # -1.474175, is negative implies cyclical behavior #R2.12.2 Decomposition of CO2 data into trend, seasonal, and #remainder sc=stl(co2, s.window=12); plot(sc) y=sc$time.series[,"trend"]; n=length(y); tim=1:n; reg1=lm(y~tim) summary(reg1) #(Intercept)=311.4468 #slope for tim =0.1092 (t stat=201.0)