• API
• FAQ
• Tools
• Archive
SHARE
TWEET # Untitled a guest Feb 16th, 2020 70 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. y=scan('D:/334/CumbriaTempMonthly.txt',skip=2,nlines=240)
2.
3. #remove the last 10% of data and plot the time series
4. y=scan('D:/334/CumbriaTempMonthly.txt',skip=2,nlines=216)
5. yts=ts(y)
6. date=seq(2000,2017+11/12,1/12)
7. plot(date,yts,main="MONTHLY TEMP IN CUMBRIA",type='p')
8.
9. #linear regression model
10. n=length(y)
11. time=seq(1:n)#^1.15
12. # where n=216 and period=12
13. month=c(rep(seq(1:12),18))
14. fmonth=as.factor(month)
15. fmonth
16.
17. ymod1=lm(y~time)
18. summary(ymod1)
19. yfit1=ymod1\$fitted
20. lines(date,yfit1,col=2)
21.
22. #Residuals of linear regression - use to check independence assumption
23. yres1=ymod1\$res
24. plot(date,yres1)
25. lines(date,yres1)
26. #are they independent? - lecture notes
27.
28. #plotting sample autocorrelation and partial autocorrelation function of residuals
29. acf(yres1)
30. pacf(yres1)
31. #numerical values
32. u=acf(yres1)
33.
34. #estimate spectral density - periodogram of residuals
36. plot(yspec\$freq,yspec\$spec,type='l',col=5)
37.
38. # harmonic regression
39. FF = 1/12 # Fundamental frequency
40. ymod2=lm(y~time+cos(2*pi*time*FF)+ sin(2*pi*time*FF))
41. summary(ymod2)
42. yfit2= ymod2\$fitted
43. plot(y, main="Temperature and harmonic regression (fundamental only)",type='p',col=2)
44. lines(yfit2,col=4)
45. lines(yfit1,col=3)
46.
47. # residuals of harmonic regression
48. ts.plot(ymod2\$res, main="Residuals after harmonic regression (fundamental only)",col=4)
49.
50. # periodogram of residuals
51. yspec=spec.pgram(ymod2\$res,pad=4,main="Periodogram of residuals after harmonic regression (fundamental only)")
52. plot(yspec\$freq,yspec\$spec,type='l',col=4,main="periodogram of residuals after harmonic regression (fundamental only)")
53.
54. # acf and pacf of harmonic residuals
55. acf(ymod2\$res,main="acf of residuals after harmonic regression (fundamental only)")
56. pacf(ymod2\$res,main="acf of residuals after harmonic regression (fundamental only)")
57.
58. # Ljung Box statistics
59. #original
60. p = rep(0, 15)
61. for (i in 1:15) {
62.   p[i] = Box.test(ymod1\$res, i, type = "Ljung-Box" )\$p.value
63. }
64. p.min = .05 #level of significance
65. plot(p,ylim=c(0,1),main="p values for Ljung-Box statistic", xlab="lag",ylab="p value")
66. abline(h=p.min,col=4,lty=2)
67. #harmonic
68. p = rep(0, 15)
69. for (i in 1:15) {
70.   p[i] = Box.test(ymod2\$res, i, type = "Ljung-Box" )\$p.value
71. }
72. p.min = .05 #level of significance
73. plot(p,ylim=c(0,1),main="p values for Ljung-Box statistic from harmonic regression", xlab="lag",ylab="p value")
74. abline(h=p.min,col=4,lty=2)
75.
76. #harmonic regression, fundamental and 2nd harmonic
77. ymod3=lm(y~time+cos(2*pi*time*FF)+ sin(2*pi*time*FF)+
78.            cos(4*pi*time*FF)+ sin(4*pi*time*FF))
79. summary(ymod3)
80. yfit3= ymod3\$fitted
81. plot(y, main="Temperature and harmonic regression - second harmonic)",type='p',col=2)
82. lines(yfit3,col=4)
83. #residuals of harmonic regression
84. ser=ts(ymod3\$res)
85. plot(ser,main="Resdiuals after harmonic fit - 2nd harmonic)",col=4)
86.
87. #periodogram of residuals
88. yspec=spec.pgram(ymod3\$res,pad=4,main="periodogram of residuals after harmonic regression - 2nd harmonic, linear scale)"
89. plot(yspec\$freq,yspec\$spec,type='l',col=4,main="periodogram of residuals after harmonic regression - with 2nd harmonic, log scale)")
90. # acf and pacf of residuals
91. acf(ymod3\$res, main="acf of residuals after harmonic regression 2nd harmonic")
92. pacf(ymod3\$res, main="acf of residuals after harmonic regression 2nd harmonic")
93. # Ljung Box statistics 2nd harmonic
94. p = rep(0, 15)
95. for (i in 1:15) {
96.   p[i] = Box.test(ymod3\$res, i, type = "Ljung-Box" )\$p.value
97. }
98. p.min = .05  # level of significance
99. plot(p, ylim=c(0,1),main="p values with Ljung-Box for residuals from harmonic regression 2nd harmonic", xlab="lag",ylab="p value")
100. abline(h=p.min,col=4,lty=2)
101.
102. #harmonic regression, 3rd harmonic
103. ymod4=lm(y~time+cos(2*pi*time*FF)+ sin(2*pi*time*FF)+
104.            cos(4*pi*time*FF)+ sin(4*pi*time*FF)+
105.            cos(6*pi*time*FF)+ sin(6*pi*time*FF))
106. summary(ymod4)
107. yfit4= ymod4\$fitted
108. plot(y, main="Cumbria Temperature and harmonic regression (with 3rd harmonic)",type='p',col=2)
109. lines(yfit4,col=4)
110. # residuals of harmonic regression
111. ser=ts(ymod4\$res)
112. plot(ser,main="Resdiuals after harmonic fit - 3rd harmonic")
113. # periodogram of residuals
114. yspec=spec.pgram(ymod4\$res,pad=4,main="periodogram of residuals after harmonic regression - 3rd harmonic, linear scale")
115. plot(yspec\$freq,yspec\$spec,type='l',col=4, main="periodogram of residuals after harmonic regression - 3rd harmonic, log scale")
116. # acf and pacf of residuals
117. acf(ymod4\$res, main="acf of residuals after harmonic regression (with 3rd harmonic)")
118. pacf(ymod4\$res, main="pacf of residuals after harmonic regression (with 3rd harmonic)")
119. # Ljung Box statistics
120. p = rep(0, 13)
121. for (i in 1:13) {
122.   p[i] = Box.test(ymod4\$res, i, type = "Ljung-Box" )\$p.value
123. }
124. p.min = .05  # level of significance
125. plot(p, ylim=c(0,1),main="p values with Ljung-Box for residuals from harmonic regression 3 harmonic", xlab="lag",ylab="p value")
126. abline(h=p.min,col=4,lty=2)
127.
128. #use periodogram to detect cycles in a time series
129. #peak in periodogram = cycle associated with a strong seasonal pattern
130. #remove trend before examining its periodogram
131. y=scan('D:/334/CumbriaTempMonthly.txt',skip=2,nlines=216)
132. n=length(y)
133. time1=seq(1:n)^1.15
134. #time1 is super linear whilst time is linear trend
135. time=seq(1:n)
136. ymod1=lm(y~time1)
137. yres1=ymod1\$res
138. ser=ts(yres1)
139. plot(ser)
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.
Top