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
  35. yspec=spec.pgram(yres1,pad=4,detrend=TRUE)
  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. OK, I Understand
Top