##################################################
#Example 2.2: matriz de dispersión (scatterplot)
mort = ts(scan("/mydata/cmort.dat"),start=1970, frequency=52)
temp = ts(scan("/mydata/temp.dat"), start=1970, frequency=52)
part = ts(scan("/mydata/part.dat"), start=1970, frequency=52)
par(mfrow=c(3,1), mar=c(3,4,3,2))
plot(mort, main="Cardiovascular Mortality", xlab="", ylab="")
plot(temp, main="Temperature", xlab="", ylab="")
plot(part, main="Particulates", xlab="", ylab="")
x11() # open another graphic device
pairs(cbind(mort, temp, part)) # scatterplot matrix
temp = temp-mean(temp)
temp2 = temp^2
trend = time(mort) # time
fit = lm(mort~ trend + temp + temp2 + part, na.action=NULL)
summary(fit) # regression results
anova(lm(mort~cbind(trend, temp, temp2, part))) # pretty ANOVA table
num = length(mort) # sample size
AIC(fit)/num - log(2*pi) # AIC (corresponds to def 2.1)
AIC(fit, k=log(num))/num - log(2*pi) # BIC (corresponds to def 2.3)
(AICc = log(sum(resid(fit)^2)/num)+(num+5)/(num-5-2)) # AICc (as in def 2.2)
#Example 2.10: suavizado por media móvil
ma5 = filter(mort, sides=2, rep(1,5)/5)
ma53 = filter(mort, sides=2, rep(1,53)/53)
plot(mort, type="p", xlab="week", ylab="mortality")
lines(ma5)
lines(ma53)
#Example 2.11: suavizado por regresión polinómica y periódica
wk = time(mort) - mean(time(mort)) # week (centered) - wk is basically t/52 ...
wk2 = wk^2 # ... because frequency=52 for mort
wk3 = wk^3
c = cos(2*pi*wk)
s = sin(2*pi*wk)
reg1 = lm(mort~wk + wk2 + wk3, na.action=NULL)
reg2 = lm(mort~wk + wk2 + wk3 + c + s, na.action=NULL)
plot(mort, type="p", ylab="mortality")
lines(fitted(reg1))
lines(fitted(reg2))
# here's the original code where mort isn't a ts object
mort = scan("/mydata/cmort.dat")
t = 1:length(mort)
t2 = t^2
t3 = t^3
c = cos(2*pi*t/52)
s = sin(2*pi*t/52)
fit1 = lm(mort~t + t2 + t3)
fit2 = lm(mort~t + t2 + t3 + c + s)
plot(t, mort)
lines(fit1$fit)
lines(fit2$fit)
#Example 2.12: suavizado kernel
plot(mort, type="p", ylab="mortality")
lines(ksmooth(time(mort), mort, "normal", bandwidth=10/52)) # or try bandwidth=5/52
lines(ksmooth(time(mort), mort, "normal", bandwidth=2))
# - original code when mort wasn't a time series with frequency=52
plot(t, mort)
lines(ksmooth(t, mort, "normal", bandwidth=10))
lines(ksmooth(t, mort, "normal", bandwidth=104))
# Red indicates a change; Fig. 2.15 shows a bandwidth of 10.
#Example 2.13: Regresión ponderada localmente y vecino más cercano
par(mfrow=c(2,1))
plot(mort,type="p", ylab="mortality", main="nearest neighbor")
lines(supsmu(time(mort), mort, span=.5))
lines(supsmu(time(mort), mort, span=.01))
plot(mort, type="p", ylab="mortality", main="lowess")
lines(lowess(mort, f=.02))
lines(lowess(mort, f=2/3))
#Example 2.14: Suavizado por splines
plot(mort, type="p", ylab="mortality")
lines(smooth.spline(time(mort), mort))
lines(smooth.spline(time(mort), mort, spar=1))
# - original code when mort wasn't a time series with frequency=52
plot(t, mort)
lines(smooth.spline(t, mort, spar=.0000001))
lines(smooth.spline(t, mort, spar=.1))
#Example 2.15: Suavizado de una serie como función de otra
par(mfrow=c(2,1))
plot(temp, mort, main="lowess")
lines(lowess(temp,mort))
plot(temp, mort, main="smoothing splines")
lines(smooth.spline(temp,mort))
Comentarios
Publicar un comentario