Advertisement
Guest User

Untitled

a guest
Apr 20th, 2014
54
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 11.19 KB | None | 0 0
  1. Date SiteSub HeatingDegreeDay MCnt MCnt_lag7 MCnt_lag9
  2. 2009-11-01 EC_BC.Z_Z 0.00 0 0 0
  3. 2009-11-02 EC_BC.Z_Z 0.00 0 0 0
  4. 2009-11-03 EC_BC.Z_Z 0.00 0 0 0
  5. 2009-11-04 EC_BC.Z_Z 0.00 0 0 0
  6. 2009-11-05 EC_BC.Z_Z 0.00 0 0 0
  7. 2009-11-06 EC_BC.Z_Z 0.00 0 0 0
  8. 2009-11-07 EC_BC.Z_Z 0.00 1 0 0
  9.  
  10. #Calculate moving average
  11. Lag0910_79 <- as.numeric(Lag0910_79$HeatingDegreeDay, Lag0910_79$MCnt7, Lag0910_79$MCnt9)
  12. Lagzoo <- as.zoo(Lag0910_79)
  13. Lagzoo_7 <- rollapply(Lagzoo, width=7, mean, na.rm=TRUE)
  14. Lagzoo_7 <- as.data.frame(Lagzoo_7)
  15.  
  16. dput(head(Lagzoo_7, 15))
  17.  
  18. structure(list(Lagzoo_7 = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
  19. 1, 1, 1, 1)), .Names = "Lagzoo_7", row.names = c("4", "5", "6",
  20. "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17",
  21. "18"), class = "data.frame")`
  22.  
  23. Lagzoo.ttr <- SMA(Lag0910_79[, "HeatingDegreeDay"], 7)
  24.  
  25. structure(list(Date = structure(c(14549, 14550, 14551, 14552,
  26. 14553, 14554, 14555, 14556, 14557, 14558, 14559, 14560, 14561,
  27. 14562, 14563, 14564, 14565, 14566, 14567, 14568, 14569, 14570,
  28. 14571, 14572, 14573, 14574, 14575, 14576, 14577, 14578, 14579,
  29. 14580, 14581, 14582, 14583, 14584, 14585, 14586, 14587, 14588,
  30. 14589, 14590, 14591, 14592, 14593, 14594, 14595, 14596, 14597,
  31. 14598, 14599, 14600, 14601, 14602, 14603, 14604, 14605, 14606,
  32. 14607, 14608, 14609, 14610, 14611, 14612, 14613, 14614, 14615,
  33. 14616, 14617, 14618, 14619, 14620, 14620, 14620, 14621, 14622,
  34. 14622, 14623, 14624, 14625, 14626, 14627, 14628, 14629, 14629,
  35. 14629, 14629, 14629, 14630, 14631, 14631, 14631, 14632, 14632,
  36. 14632, 14632, 14632, 14632, 14632, 14633, 14633, 14633, 14634,
  37. 14634, 14634, 14634, 14635, 14635, 14635, 14635, 14636, 14636,
  38. 14636, 14636, 14636, 14636, 14637, 14637, 14637, 14638, 14638,
  39. 14638, 14639, 14639, 14640, 14641, 14642, 14643, 14643, 14644,
  40. 14645, 14646, 14647, 14648, 14649, 14650, 14651, 14652, 14653,
  41. 14654, 14655, 14656, 14657, 14658, 14659, 14660, 14661, 14661,
  42. 14662, 14663, 14663, 14664, 14665, 14666, 14667, 14668, 14669,
  43. 14669, 14670, 14671, 14672, 14673, 14674, 14675, 14675, 14676,
  44. 14677, 14678, 14678, 14679, 14680, 14681, 14681, 14681, 14682,
  45. 14682, 14682, 14683, 14684, 14685, 14686, 14687, 14688, 14689,
  46. 14689, 14690, 14691, 14692, 14693, 14694, 14694, 14694, 14695,
  47. 14696, 14697, 14698, 14699, 14700, 14701, 14702, 14703, 14703,
  48. 14703, 14703, 14704, 14704, 14705, 14706), class = "Date"), SiteSub = structure(c(1L,
  49. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  50. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  51. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  52. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  53. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  54. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  55. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  56. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  57. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  58. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  59. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  60. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  61. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = "EC_BC.Z_Z", class = "factor"),
  62. HeatingDegreeDay = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L,
  63. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  64. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  65. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  66. 1L, 1L, 1L, 1L, 1L, 1L, 1L, NA, 1L, 1L, 3L, 4L, 9L, 11L,
  67. 14L, 15L, 12L, 13L, 17L, 16L, 16L, 16L, 10L, 8L, 8L, 7L,
  68. 6L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  69. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  70. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  71. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  72. 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  73. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 5L, 1L, 1L, 1L, 1L, 1L, 1L,
  74. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  75. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  76. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("0.00",
  77. "0.02", "0.05", "0.14", "0.32", "0.50", "0.89", "0.96", "0.98",
  78. "1.02", "1.04", "1.30", "1.40", "1.49", "1.50", "1.58", "1.86"
  79. ), class = "factor"), MCnt = structure(c(1L, 1L, 1L, 1L,
  80. 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  81. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  82. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  83. 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  84. 1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 1L,
  85. 2L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
  86. 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
  87. 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
  88. 1L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L,
  89. 1L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 1L,
  90. 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
  91. 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L,
  92. 2L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 1L,
  93. 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("0", "1"), class = "factor"),
  94. MCnt_lag7 = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  95. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  96. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  97. 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  98. 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 2L,
  99. 2L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
  100. 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
  101. 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
  102. 2L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L,
  103. 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 2L,
  104. 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L,
  105. 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
  106. 1L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
  107. 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, NA, NA, NA,
  108. NA, NA, NA, NA), .Label = c("0", "1"), class = "factor"),
  109. MCnt_lag9 = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  110. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  111. 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  112. 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  113. 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 2L,
  114. 1L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
  115. 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
  116. 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
  117. 2L, 1L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L,
  118. 1L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 1L,
  119. 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
  120. 2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L,
  121. 2L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L,
  122. 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, NA, NA, NA, NA, NA,
  123. NA, NA, NA, NA), .Label = c("0", "1"), class = "factor")), .Names = c("Date",
  124. "SiteSub", "HeatingDegreeDay", "MCnt", "MCnt_lag7", "MCnt_lag9"
  125. ), row.names = c("1", "2", "3", "4", "5", "6", "7", "8", "9",
  126. "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "20",
  127. "21", "22", "23", "24", "25", "26", "27", "28", "29", "30", "31",
  128. "32", "33", "34", "35", "36", "37", "38", "39", "40", "41", "42",
  129. "43", "44", "45", "46", "47", "48", "49", "50", "51", "52", "53",
  130. "54", "55", "56", "57", "58", "59", "60", "61", "62", "63", "64",
  131. "65", "66", "67", "68", "69", "70", "71", "72", "73", "74", "75",
  132. "76", "77", "78", "79", "80", "81", "82", "83", "84", "85", "86",
  133. "87", "88", "89", "90", "91", "92", "93", "94", "95", "96", "97",
  134. "98", "99", "100", "101", "102", "103", "104", "105", "106",
  135. "107", "108", "109", "110", "111", "112", "113", "114", "115",
  136. "116", "117", "118", "119", "120", "121", "122", "123", "124",
  137. "125", "126", "127", "128", "129", "130", "131", "132", "133",
  138. "134", "135", "136", "137", "138", "139", "140", "141", "142",
  139. "143", "144", "145", "146", "147", "148", "149", "150", "151",
  140. "152", "153", "154", "155", "156", "157", "158", "159", "160",
  141. "161", "162", "163", "164", "165", "166", "167", "168", "169",
  142. "170", "171", "172", "173", "174", "175", "176", "177", "178",
  143. "179", "180", "181", "182", "183", "184", "185", "186", "187",
  144. "188", "189", "190", "191", "192", "193", "194", "195", "196",
  145. "197", "198", "199", "200", "201", "202", "203", "204", "205",
  146. "206", "207", "208"), class = "data.frame")
  147.  
  148. library(zoo)
  149.  
  150. DF2 <- DF
  151. DF2[3:6] <- lapply(DF2[3:6], function(x) as.numeric(as.character(x)))
  152. m <- as.matrix(DF2[3:6])
  153. rmean <- rollapplyr(m, 7, mean, na.rm = TRUE, fill = NA) # mean matrix
  154.  
  155. corNA <- function(x) {
  156. x <- na.omit(x[, 1:2])
  157. if (nrow(x) < 2 || sd(x[,1]) == 0 || sd(x[,2]) == 0) return(NA)
  158. cor(x[, 1], x[,2])
  159. }
  160.  
  161. rcor <- rollapplyr(m, 7, corNA, by.column = FALSE, fill = NA) # vector of cors
  162.  
  163. DF3 <- data.frame(DF2, rmean, rcor) # put it all together
  164.  
  165. z <- read.zoo(DF2[-2], aggregate = mean) # can omit aggregate=mean if dates are unique
  166.  
  167. zmean <- rollapplyr(z, 7, mean, na.rm = TRUE, fill = NA) # means
  168. zcor <- rollapplyr(z, 7, corNA, by.column = FALSE, fill = NA) # cors
  169.  
  170. z2 <- merge(z, zmean, zcor) # omit this if separate objects are ok
  171.  
  172. # Convert dates to days
  173. aa = as.Date(x$Date)
  174. x$Date = as.numeric(aa - aa[1])
  175.  
  176. # I think it's easier to get rid of the factors
  177. factor2number = function(x) as.numeric(as.character(x))
  178. x[,3:6] = apply(x[,3:6],2,factor2number)
  179.  
  180. # A rolling mean function
  181. rollmean_r = function(x,y,width) {
  182. out = numeric(length(x))
  183. for( i in seq_along(x) ) {
  184. window = x >= (x[i]-width) & x <= (x[i]+width)
  185. out[i] = .Internal(mean( y[window] ))
  186. }
  187. return(out)
  188. }
  189.  
  190. # Calculate the rolling means
  191. x[,3:6] = apply(x[3:6], 2, function(y) rollmean_r(x$Date,y,7) )
  192. x
  193. # Date SiteSub HeatingDegreeDay MCnt MCnt_lag7 MCnt_lag9
  194. # 1 0 EC_BC.Z_Z 0 0.12500000 0 0
  195. # 2 1 EC_BC.Z_Z 0 0.11111111 0 0
  196. # 3 2 EC_BC.Z_Z 0 0.10000000 0 0
  197. # 4 3 EC_BC.Z_Z 0 0.09090909 0 0
  198. # 5 4 EC_BC.Z_Z 0 0.08333333 0 0
  199. # 6 5 EC_BC.Z_Z 0 0.07692308 0 0
  200. # 7 6 EC_BC.Z_Z 0 0.07142857 0 0
  201. # 8 7 EC_BC.Z_Z 0 0.06666667 0 0
  202. # 9 8 EC_BC.Z_Z 0 0.07142857 0 0
  203. # 10 9 EC_BC.Z_Z 0 0.07692308 0 0
  204. # 11 10 EC_BC.Z_Z 0 0.08333333 0 0
  205. # 12 11 EC_BC.Z_Z 0 0.09090909 0 0
  206. # 13 12 EC_BC.Z_Z 0 0.10000000 0 0
  207. # 14 13 EC_BC.Z_Z 0 0.11111111 0 0
  208. # 15 14 EC_BC.Z_Z 0 0.00000000 0 0
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement