MATD51-Aula 07: Manchas Solares como soma de cossenos
Kim Samejima
3/28/2021
Análise da sazonalidade da série de manchas solares.
Carregando os dados da série de manchas solares:
db_sunspots=read.table("../Dados/manchas.csv", sep=",",h=T)
T=length(db_sunspots$manchas)
t=1:TGráfico da série: plot, acf e pacf
ppacf(db_sunspots$manchas, new.window = FALSE, main="Série de manchas solares",
evaldist = TRUE)
Espectro da Série
spec_sunspots=spec0(db_sunspots$manchas, id.type = "freq", nfreqs = 4, new.window = FALSE)
Identificando a periodicidade da série original:
w1=spec_sunspots$spectrum$freq[spec_sunspots$pos.maxspecs[1]] # maior energia "primeiro pico do spectro"
w2=spec_sunspots$spectrum$freq[spec_sunspots$pos.maxspecs[2]] # 2a maior energia "segundo pico do spectro"
w3=spec_sunspots$spectrum$freq[spec_sunspots$pos.maxspecs[3]] # 3a maior energia "terceiro pico do spectro"
w4=spec_sunspots$spectrum$freq[spec_sunspots$pos.maxspecs[4]] # 4a maior energia "quarto pico do spectro" A partir do perÃodograma, vemos que as quarto maiores periodicidades são: 11.25 (frequência 0.0888889), 90 (frequência 0.0111111), 12 (frequência 0.0833333) e 10 (frequência 0.1).
Ajustando o modelo periódico para a maior periodicidade
#w1=0.088888889
# Construindo as variaveis z1 e z2
z1 = cos(2*pi*w1*t)
z2 = sin(2*pi*w1*t)
# Ajustando o modelo para encontrar beta1 e beta2
fit_cosw1 <- lm(db_sunspots$manchas~z1+z2)
summary(fit_cosw1)##
## Call:
## lm(formula = db_sunspots$manchas ~ z1 + z2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -52.738 -18.296 -5.362 12.472 118.418
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 45.158 2.196 20.561 < 2e-16 ***
## z1 26.730 3.104 8.611 4.34e-15 ***
## z2 -4.837 3.107 -1.557 0.121
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 29.13 on 173 degrees of freedom
## Multiple R-squared: 0.3062, Adjusted R-squared: 0.2982
## F-statistic: 38.17 on 2 and 173 DF, p-value: 1.851e-14# beta1 e beta2
b1b2=fit_cosw1$coefficients
# encontrando phi1 e A1
phi1=atan(-b1b2[2]/b1b2[1])
A1=b1b2[1]/cos(phi1)
print(paste("w1=",round(w1,2),";","A=",round(A1,2),";","Phi=",round(phi1,2)),digits = 2)## [1] "w1= 0.09 ; A= 52.48 ; Phi= -0.53"Gráfico da série e do resÃduo
par(mfrow=c(1,2),cex.axis=0.8, mar=c(2,2,4,2))
plot.ts(db_sunspots$manchas, col="darkgrey", lwd=2,ylab="",
main="Série de manchas solares \n e modelo sazonal ajustado")
lines(fitted(fit_cosw1), col="red",lty=2)
resid_sunspots=db_sunspots$manchas-fitted(fit_cosw1);
plot.ts(resid_sunspots,ylab="",main="ResÃduos do ajuste sazonal")
Com este modelo, terÃamos que a série de manchas solares (\(X_t\)) poderia ser modelada pelo modelo: \[ X_t = A_1cos(2\pi\omega_1 t +\phi_1)+\epsilon_{1,t}, \] com \(\widehat{A_1}\)=52.4756216 e \(\widehat{\phi_1}\)=-0.5344579.
Ajustando o modelo periódico para a \(2^a\) maior periodicidade
Podemos, iterativamente, modelar o resÃduo do modelo anterior para captar uma segunda periodicidade. Vejamos como fica o perÃodograma do resÃduo \(\epsilon_{1,t}\):
# perÃodograma do resÃduo
spec_resid1=spec0(resid_sunspots, id.type = "period", nfreqs = 4, new.window = FALSE)
w2=spec_resid1$spectrum$freq[spec_resid1$pos.maxspecs[1]] # 2a maior energia "primeiro pico do spectro do primeiro resÃduo" Assim, para modelar a periodicidade do resÃduo, farÃamos o ajuste de forma análoga ao modelo anterior:
# Construindo as variaveis z3 e z4
z3 = cos(2*pi*w2*t)
z4 = sin(2*pi*w2*t)
# Ajustando o modelo para encontrar beta3 e beta4
fit_cosw2 <- lm(resid_sunspots~z3+z4)
summary(fit_cosw2)##
## Call:
## lm(formula = resid_sunspots ~ z3 + z4)
##
## Residuals:
## Min 1Q Median 3Q Max
## -48.209 -15.831 -4.897 9.601 109.628
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.01517 2.01333 0.008 0.994
## z3 2.35264 2.87908 0.817 0.415
## z4 15.97931 2.81508 5.676 5.7e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 26.7 on 173 degrees of freedom
## Multiple R-squared: 0.1599, Adjusted R-squared: 0.1502
## F-statistic: 16.47 on 2 and 173 DF, p-value: 2.842e-07# beta3 e beta4
b3b4=fit_cosw2$coefficients
# encontrando phi2 e A2
phi2=atan(-b3b4[2]/b3b4[1])
A2=b3b4[1]/cos(phi2)
print(paste("w2=",round(w2,2),";","A2=",round(A2,2),";","Phi2=",round(phi2,2)),digits = 2)## [1] "w2= 0.01 ; A2= 2.35 ; Phi2= -1.56"Gráfico da série e do resÃduo
par(mfrow=c(1,2),cex.axis=0.8, mar=c(2,2,4,2))
plot.ts(resid_sunspots, col="darkgrey", lwd=2,ylab="",
main="ResÃduo 1 \n e modelo sazonal ajustado")
lines(fitted(fit_cosw2), col="red",lty=2)
resid_sunspots2=resid_sunspots-fitted(fit_cosw2);
plot.ts(resid_sunspots2,ylab="",main="ResÃduos do segundo ajuste sazonal")
A série \(\epsilon_{1,t}\) pode portanto ser escrita como: \[
\epsilon_{1,t} = A_2cos(2\pi\omega_2 t +\phi_2)+\epsilon_{2,t},
\] com \(\widehat{A_2}\)=2.3526928 e \(\widehat{\phi_2}\)=-1.5643489.
Ou, se quisermos olhara para a série original, então a série de manchas solares (\(X_t\)) pode ser modelada como: \[
X_t = A_1cos(2\pi\omega_1 t +\phi_1)+A_2cos(2\pi\omega_2 t +\phi_2)+\epsilon_{2,t}.
\]
plot.ts(db_sunspots$manchas, col="darkgrey", lwd=2,ylab="",
main="Série de manchas solares \n e modelo sazonal ajustado com duas componentes")
lines(fitted(fit_cosw1)+fitted(fit_cosw2), col="red",lty=2,lwd=2)
Ajustando o modelo periódico para a \(3^a\) maior periodicidade
# perÃodograma do resÃduo2
spec_resid2=spec0(resid_sunspots2, id.type = "period", nfreqs = 4, new.window = FALSE)# kim
w3=spec_resid2$spectrum$freq[spec_resid2$pos.maxspecs[1]] # 3a maior energia "primeiro pico do spectro do segundo resÃduo"
#w3=0.1
# Construindo as variaveis z5 e z6
z5 = cos(2*pi*w3*t)
z6 = sin(2*pi*w3*t)
# Ajustando o modelo para encontrar beta5 e beta6
fit_cosw3 <- lm(resid_sunspots2~z5+z6)
summary(fit_cosw3)##
## Call:
## lm(formula = resid_sunspots2 ~ z5 + z6)
##
## Residuals:
## Min 1Q Median 3Q Max
## -46.186 -13.292 -2.446 8.200 103.756
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.02203 1.84063 0.012 0.99047
## z5 12.92051 2.60028 4.969 1.61e-06 ***
## z6 7.83020 2.60510 3.006 0.00304 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 24.41 on 173 degrees of freedom
## Multiple R-squared: 0.1638, Adjusted R-squared: 0.1542
## F-statistic: 16.95 on 2 and 173 DF, p-value: 1.899e-07# beta5 e beta6
b5b6=fit_cosw3$coefficients
# encontrando phi3 e A3
phi3=atan(-b5b6[2]/b5b6[1])
A3=b5b6[1]/cos(phi3)
print(paste("w3=",round(w3,2),";","A3=",round(A3,2),";","Phi3=",round(phi3,2)),digits = 2)## [1] "w3= 0.1 ; A3= 12.92 ; Phi3= -1.57"Gráfico da série e do resÃduo
par(mfrow=c(1,2),cex.axis=0.8, mar=c(2,2,4,2))
plot.ts(resid_sunspots2, col="darkgrey", lwd=2,ylab="",
main="ResÃduo 2 \n e modelo sazonal ajustado")
lines(fitted(fit_cosw3), col="red",lty=2)
resid_sunspots3=resid_sunspots2-fitted(fit_cosw3);
plot.ts(resid_sunspots3,ylab="",main="ResÃduos do segundo ajuste sazonal")
plot.ts(db_sunspots$manchas, col="darkgrey", lwd=2,ylab="",
main="Série de manchas solares \n e modelo sazonal ajustado com três componentes")
lines(fitted(fit_cosw1)+fitted(fit_cosw2)+fitted(fit_cosw3), col="red",lty=2,lwd=2)
A série \(\epsilon_{2,t}\) pode portanto ser escrita como: \[
\epsilon_{2,t} = A_3cos(2\pi\omega_3 t +\phi_3)+\epsilon_{3,t},
\] com \(\widehat{A_3}\)=12.9205268 e \(\widehat{\phi_3}\)=-1.5690914.
Ou, se quisermos olhara para a série original, então a série de manchas solares (\(X_t\)) pode ser modelada como: \[
X_t = A_1cos(2\pi\omega_1 t +\phi_1)+A_2cos(2\pi\omega_2 t +\phi_2)+A_3cos(2\pi\omega_3 t +\phi_3)+\epsilon_{3,t}.
\]
Ajustando o modelo periódico para a \(4^a\) maior periodicidade
# perÃodograma do resÃduo4
spec_resid3=spec0(resid_sunspots3, id.type = "period", nfreqs = 4, new.window = FALSE)# kim
w4=spec_resid3$spectrum$freq[spec_resid3$pos.maxspecs[1]] # 4a maior energia "primeiro pico do spectro do terceiro resÃduo"
#w4=0.08333333
# Construindo as variaveis z7 e z8
z7 = cos(2*pi*w4*t)
z8 = sin(2*pi*w4*t)
# Ajustando o modelo para encontrar beta7 e beta8
fit_cosw4 <- lm(resid_sunspots3~z7+z8)
summary(fit_cosw4)##
## Call:
## lm(formula = resid_sunspots3 ~ z7 + z8)
##
## Residuals:
## Min 1Q Median 3Q Max
## -39.321 -13.585 -1.322 8.927 92.103
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.2402 1.6971 -0.142 0.8876
## z7 -5.9773 2.3998 -2.491 0.0137 *
## z8 11.8935 2.3998 4.956 1.7e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 22.51 on 173 degrees of freedom
## Multiple R-squared: 0.15, Adjusted R-squared: 0.1401
## F-statistic: 15.26 on 2 and 173 DF, p-value: 7.887e-07# beta7 e beta8
b7b8=fit_cosw4$coefficients
# encontrando phi4 e A4
phi4=atan(-b7b8[2]/b7b8[1])
A4=b7b8[1]/cos(phi4)
print(paste("w4=",round(w4,2),";","A4=",round(A4,2),";","Phi4=",round(phi4,2)),digits = 2)## [1] "w4= 0.08 ; A4= -5.98 ; Phi4= -1.53"Gráfico da série e do resÃduo
par(mfrow=c(1,2),cex.axis=0.8, mar=c(2,2,4,2))
plot.ts(resid_sunspots3, col="darkgrey", lwd=2,ylab="",
main="ResÃduo 2 \n e modelo sazonal ajustado")
lines(fitted(fit_cosw4), col="red",lty=2)
resid_sunspots4=resid_sunspots3-fitted(fit_cosw4);
plot.ts(resid_sunspots4,ylab="",main="ResÃduos do segundo ajuste sazonal")
plot.ts(db_sunspots$manchas, col="darkgrey", lwd=2,ylab="",
main="Série de manchas solares \n e modelo sazonal ajustado com quatro componentes")
lines(fitted(fit_cosw1)+fitted(fit_cosw2)+fitted(fit_cosw3)+fitted(fit_cosw4), col="red",lty=2,lwd=2)
# perÃodograma do resÃduo4
spec_resid4=spec0(resid_sunspots4, id.type = "period", nfreqs = 4, new.window = FALSE)
A série \(\epsilon_{3,t}\) pode portanto ser escrita como: \[
\epsilon_{3,t} = A_4cos(2\pi\omega_4 t +\phi_4)+\epsilon_{4,t},
\] com \(\widehat{A_4}\)=-5.9821391 e \(\widehat{\phi_4}\)=-1.5306255.
…e assim sucessivamente
par(mfrow=c(2,2))
plot.ts(db_sunspots$manchas, col="darkgrey", lwd=2,ylab="",
main="Série de manchas solares \n e modelo sazonal ajustado")
lines(fitted(fit_cosw1), col="orange",lty=1,lwd=2)
plot.ts(db_sunspots$manchas, col="darkgrey", lwd=2,ylab="",
main="Série de manchas solares \n e modelo sazonal ajustado")
lines(fitted(fit_cosw1)+fitted(fit_cosw2), col="red",lty=1,lwd=2)
plot.ts(db_sunspots$manchas, col="darkgrey", lwd=2,ylab="",
main="Série de manchas solares \n e modelo sazonal ajustado")
lines(fitted(fit_cosw1)+fitted(fit_cosw2)+fitted(fit_cosw3), col="blue",lty=1,lwd=2)
plot.ts(db_sunspots$manchas, col="darkgrey", lwd=2,ylab="",
main="Série de manchas solares \n e modelo sazonal ajustado")
lines(fitted(fit_cosw1)+fitted(fit_cosw2)+fitted(fit_cosw3)+fitted(fit_cosw4), col="purple",lty=1,lwd=2)
par(mfrow=c(2,2))
plot.ts(fitted(fit_cosw1), col="orange",lty=1,lwd=2, ylab="", main="Ajuste 1")
plot.ts(fitted(fit_cosw1)+fitted(fit_cosw2), col="red",lty=1,lwd=2, ylab="", main="Ajuste 1+2")
plot.ts(fitted(fit_cosw1)+fitted(fit_cosw2)+fitted(fit_cosw3), col="blue",lty=1,lwd=2, ylab="", main="Ajuste 1+2+3")
plot.ts(fitted(fit_cosw1)+fitted(fit_cosw2)+fitted(fit_cosw3)+fitted(fit_cosw4), col="purple",lty=1,lwd=2, ylab="", main="Ajuste 1+2+3+4")
A série de manchas solares (\(X_t\)) pode ser modelada como: \[ X_t = A_1cos(2\pi\omega_1 t +\phi_1)+A_2cos(2\pi\omega_2 t +\phi_2)+A_3cos(2\pi\omega_3 t +\phi_3)+A_4cos(2\pi\omega_4 t +\phi_4)+\epsilon_{4,t}, \] com
com as periodididades (11.25, 90, 10, 12) modeladas pelos parâmetros \((\widehat{A_1},\widehat{A_2},\widehat{A_3},\widehat{A_4})\)=(52.48, 2.35, 12.92, -5.98) e \((\widehat{\phi_1},\widehat{\phi_2},\widehat{\phi_3},\widehat{\phi_4})\)=(-0.53, -1.56, -1.57, -1.53).
Modelando de forma simultânea
Os ajustes dos coeficientes do modelo acima foram feitos de forma independente, isto é, cada modelo subsequente foi ajustado com o resÃduo do modelo anterior. O ideal é ajustar o modelo de forma simultânea, para se encontrar o ótimo global para todos os parâmetros estimados.
Identificando em quais frequências ocorrem as maiores energias:
nfreqs=12
spec_sunspots=spec0(db_sunspots$manchas, id.type = "period", nfreqs = nfreqs, new.window = FALSE)
w=spec_sunspots$spectrum$freq[spec_sunspots$pos.maxspecs[1:nfreqs]]
z1=cos(2*pi*cbind(t)%*%rbind(w))
z2=sin(2*pi*cbind(t)%*%rbind(w))Ajustando o modelo periódico com diversos coeficientes simultaneamente
# Construindo as variaveis z1 e z2
# Ajustando o modelo para encontrar beta1 ... beta12
fit_6_w <- lm(db_sunspots$manchas~z1+z2)
summary(fit_6_w)##
## Call:
## lm(formula = db_sunspots$manchas ~ z1 + z2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -34.143 -9.220 -0.093 9.387 55.304
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 44.9569 1.1641 38.618 < 2e-16 ***
## z11 26.5691 1.6421 16.180 < 2e-16 ***
## z12 1.9253 1.6748 1.150 0.252131
## z13 -6.1287 1.6439 -3.728 0.000272 ***
## z14 12.5429 1.6396 7.650 2.20e-12 ***
## z15 -9.5396 1.6736 -5.700 6.10e-08 ***
## z16 -9.4778 1.6382 -5.785 4.03e-08 ***
## z17 2.5355 1.6383 1.548 0.123804
## z18 0.8448 1.6407 0.515 0.607363
## z19 -1.2563 1.6388 -0.767 0.444529
## z110 5.0727 1.6720 3.034 0.002842 **
## z111 0.5093 1.6509 0.309 0.758115
## z112 -1.4284 1.6460 -0.868 0.386867
## z21 -4.5411 1.6451 -2.760 0.006488 **
## z22 16.0308 1.6173 9.912 < 2e-16 ***
## z23 11.7254 1.6439 7.133 3.82e-11 ***
## z24 8.0109 1.6463 4.866 2.84e-06 ***
## z25 -4.7891 1.6184 -2.959 0.003584 **
## z26 -2.9325 1.6457 -1.782 0.076773 .
## z27 -10.0237 1.6450 -6.094 8.79e-09 ***
## z28 -7.3329 1.6459 -4.455 1.62e-05 ***
## z29 7.7788 1.6464 4.725 5.23e-06 ***
## z210 -6.2987 1.6199 -3.888 0.000151 ***
## z211 -5.4631 1.6385 -3.334 0.001076 **
## z212 7.0274 1.6424 4.279 3.32e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.33 on 151 degrees of freedom
## Multiple R-squared: 0.8321, Adjusted R-squared: 0.8055
## F-statistic: 31.19 on 24 and 151 DF, p-value: < 2.2e-16# beta1 ... beta12
b=fit_6_w$coefficients[-1]
# encontrando phi e A
phi=A=numeric()
for( i in 1:nfreqs){
phi[i]=atan(-b[nfreqs+i]/b[i])
A[i]=b[i]/cos(phi[i])
}Periodididades (11.25, 90, 12, 10, 60, 8.57, 8.18, 10.59, 9.47, 45, 15, 12.86) modeladas pelos parâmetros \((\widehat{A_1},\dots,\widehat{A_{12}})\)=(26.95, 16.15, -13.23, 14.88, -10.67, -9.92, 10.34, 7.38, -7.88, 8.09, 5.49, -7.17) e \((\widehat{\phi_1},\dots,\widehat{\phi_{12}})\)=(0.17, -1.45, 1.09, -0.57, -0.47, -0.3, 1.32, 1.46, 1.41, 0.89, 1.48, 1.37).
Gráfico da série e do resÃduo
par(mfrow=c(2,2),cex.axis=0.8, mar=c(2,2,4,2))
plot.ts(db_sunspots$manchas, col="darkblue", lwd=2,ylab="",
main="Série de manchas solares")#, ylim=c(0,150))
plot.ts(fitted(fit_6_w), col="orange", lwd=2,ylab="",
main="Modelo sazonal ajustado")#, ylim=c(0,150))
plot.ts(db_sunspots$manchas, col="darkgrey", lwd=2,ylab="",
main="Série de manchas solares \n e modelo sazonal ajustado")
lines(fitted(fit_6_w), col="red",lty=2, lwd=2)
resid_sunspots=db_sunspots$manchas-fitted(fit_6_w);
plot.ts(resid_sunspots,ylab="",main="ResÃduos do ajuste sazonal")