Statistics and Data Visualizations in Climate Science
with R and Python
A Cambridge University Press Book by
SSP Shen and GR North
Version 1.0 released in July 2023 and coded by Dr. Samuel Shen,
Distinguished Professor
San Diego State University, California, USA
Email: sshen@sdsu.edu
Chapter 7: Introduction to Time Series
R Code for Fig. 7.1: Mauna Loa CO2 March 1958-July 2020
#setwd("/Users/sshen/climstats")
co2m = read.table("data/co2m.txt", header = TRUE)
dim(co2m)## [1] 749 7# [1] 749 7
co2m[1:3,]## year mon date average interpolated trend days
## 1 1958 3 1958.208 315.71 315.71 314.62 -1
## 2 1958 4 1958.292 317.45 317.45 315.29 -1
## 3 1958 5 1958.375 317.50 317.50 314.71 -1#year mon date average interpolated trend days
#1 1958 3 1958.208 315.71 315.71 314.62 -1
#2 1958 4 1958.292 317.45 317.45 315.29 -1
#3 1958 5 1958.375 317.50 317.50 314.71 -1
mon = co2m[,3]
co2 = co2m[,5]
#setEPS() # save the .eps figure file to the working directory
#postscript("fig0701.eps", height = 8, width = 10)
par(mar=c(2.2,4.5,2,0.5))
plot(mon, co2, type="l", col="red",
main =
expression(paste("Monthly Atmospheric ",
CO[2]," at Mauna Loa Observatory")),
xlab ="",
ylab="Parts Per Million [ppm]",
cex.axis =1.4, cex.lab=1.4, cex.main= 1.5)
text(1985, 410, "Scripps Institution of Oceanography",
cex=1.5)
text(1985, 400, "NOAA Global Monitoring Laboratory",
cex=1.5)
lines(mon, co2m[,6]) #plot the trend data
#dev.off()
R Plot Fig. 7.2: Time Series Decomposition
co2.ts = ts(co2, start=c(1958,3), end=c(2020,7),
frequency =12)
co2.ts #Display time series with month and year## Jan Feb Mar Apr May Jun Jul Aug Sep Oct
## 1958 315.71 317.45 317.50 317.10 315.86 314.93 313.20 312.66
## 1959 315.62 316.38 316.71 317.72 318.29 318.15 316.54 314.80 313.84 313.26
## 1960 316.43 316.97 317.58 319.02 320.03 319.59 318.18 315.91 314.16 313.83
## 1961 316.93 317.70 318.54 319.48 320.58 319.77 318.57 316.79 314.80 315.38
## 1962 317.94 318.56 319.68 320.63 321.01 320.55 319.58 317.40 316.26 315.42
## 1963 318.74 319.08 319.86 321.39 322.25 321.47 319.74 317.77 316.21 315.99
## 1964 319.57 320.07 320.73 321.77 322.25 321.89 320.44 318.70 316.70 316.79
## 1965 319.44 320.44 320.89 322.13 322.16 321.87 321.39 318.81 317.81 317.30
## 1966 320.62 321.59 322.39 323.87 324.01 323.75 322.39 320.37 318.64 318.10
## 1967 322.07 322.50 323.04 324.42 325.00 324.09 322.55 320.92 319.31 319.31
## 1968 322.57 323.15 323.89 325.02 325.57 325.36 324.14 322.03 320.41 320.25
## 1969 324.00 324.42 325.64 326.66 327.34 326.76 325.88 323.67 322.38 321.78
## 1970 325.03 325.99 326.87 328.13 328.07 327.66 326.35 324.69 323.10 323.16
## 1971 326.17 326.68 327.18 327.78 328.92 328.57 327.34 325.46 323.36 323.57
## 1972 326.77 327.63 327.75 329.72 330.07 329.09 328.05 326.32 324.93 325.06
## 1973 328.54 329.56 330.30 331.50 332.48 332.07 330.87 329.31 327.51 327.18
## 1974 329.35 330.71 331.48 332.65 333.19 332.16 331.07 329.12 327.32 327.28
## 1975 330.73 331.46 331.90 333.17 333.94 333.45 331.97 329.95 328.50 328.34
## 1976 331.59 332.75 333.52 334.64 334.77 334.00 333.06 330.68 328.95 328.75
## 1977 332.66 333.13 334.95 336.13 336.93 336.17 334.88 332.56 331.29 331.27
## 1978 334.95 335.25 336.66 337.69 338.03 338.01 336.41 334.41 332.37 332.41
## 1979 336.14 336.69 338.27 338.95 339.21 339.26 337.54 335.75 333.98 334.19
## 1980 337.90 338.34 340.01 340.93 341.48 341.33 339.40 337.70 336.19 336.15
## 1981 339.29 340.55 341.61 342.53 343.03 342.54 340.78 338.44 336.95 337.08
## 1982 340.96 341.73 342.81 343.97 344.63 343.79 342.32 340.09 338.28 338.29
## 1983 341.68 342.90 343.33 345.25 346.03 345.63 344.19 342.27 340.35 340.38
## 1984 344.10 344.79 345.52 346.84 347.63 346.98 345.53 343.55 341.40 341.67
## 1985 345.21 346.16 347.74 348.34 349.06 348.38 346.71 345.02 343.27 343.13
## 1986 346.56 347.28 348.01 349.77 350.38 349.93 348.16 346.08 345.22 344.51
## 1987 348.52 348.73 349.73 351.31 352.09 351.53 350.11 348.08 346.52 346.59
## 1988 350.39 351.64 352.40 353.69 354.21 353.72 352.69 350.40 348.92 349.13
## 1989 352.91 353.27 353.96 355.64 355.86 355.37 353.99 351.81 350.05 350.25
## 1990 353.80 355.04 355.73 356.32 357.32 356.34 354.84 353.01 351.31 351.62
## 1991 354.84 355.73 357.23 358.66 359.13 358.13 356.19 353.85 352.25 352.35
## 1992 356.25 357.11 357.86 359.09 359.59 359.33 357.01 354.94 352.95 353.32
## 1993 357.00 357.31 358.47 359.27 360.19 359.52 357.33 355.64 354.03 354.12
## 1994 358.24 358.92 359.99 361.23 361.65 360.81 359.38 357.46 355.73 356.08
## 1995 359.92 360.86 361.83 363.30 363.69 363.19 361.64 359.12 358.17 357.99
## 1996 362.07 363.24 364.17 364.57 365.13 364.92 363.55 361.38 359.54 359.58
## 1997 363.09 364.03 364.51 366.35 366.64 365.59 364.31 362.25 360.29 360.82
## 1998 365.27 365.98 367.24 368.66 369.42 368.99 367.82 365.95 364.02 364.40
## 1999 368.18 369.07 369.68 370.99 370.96 370.30 369.45 366.90 364.81 365.37
## 2000 369.29 369.55 370.60 371.82 371.58 371.70 369.86 368.13 367.00 367.03
## 2001 370.59 371.51 372.43 373.37 373.85 373.22 371.50 369.61 368.18 368.45
## 2002 372.53 373.20 374.12 375.02 375.76 375.52 374.01 371.85 370.75 370.55
## 2003 374.88 375.64 376.45 377.73 378.60 378.28 376.70 374.38 373.17 373.15
## 2004 377.00 377.87 378.88 380.35 380.62 379.69 377.47 376.01 374.25 374.46
## 2005 378.46 379.73 380.77 382.29 382.45 382.21 380.74 378.74 376.70 377.00
## 2006 381.38 382.20 382.67 384.61 385.03 384.05 382.46 380.41 378.85 379.13
## 2007 382.89 383.90 384.58 386.50 386.56 386.10 384.50 381.99 380.96 381.12
## 2008 385.52 385.82 386.03 387.21 388.54 387.76 386.37 384.09 383.18 382.99
## 2009 386.94 387.48 388.82 389.55 390.14 389.48 388.03 386.11 384.74 384.43
## 2010 388.71 390.20 391.17 392.46 393.00 392.15 390.20 388.35 386.85 387.24
## 2011 391.33 391.86 392.60 393.25 394.19 393.74 392.51 390.13 389.08 389.00
## 2012 393.12 393.86 394.40 396.18 396.74 395.71 394.36 392.39 391.11 391.05
## 2013 395.55 396.80 397.43 398.41 399.78 398.60 397.32 395.20 393.45 393.70
## 2014 397.85 398.01 399.77 401.38 401.78 401.25 399.10 397.03 395.38 396.03
## 2015 399.98 400.28 401.54 403.28 403.96 402.80 401.31 398.93 397.63 398.29
## 2016 402.56 404.12 404.87 407.45 407.72 406.83 404.41 402.27 401.05 401.59
## 2017 406.17 406.46 407.22 409.04 409.69 408.88 407.12 405.13 403.37 403.63
## 2018 407.96 408.32 409.41 410.24 411.24 410.79 408.71 406.99 405.51 406.00
## 2019 410.83 411.75 411.97 413.33 414.64 413.93 411.74 409.95 408.54 408.52
## 2020 413.39 414.11 414.51 416.21 417.07 416.38 414.38
## Nov Dec
## 1958 313.33 314.67
## 1959 314.80 315.58
## 1960 315.00 316.19
## 1961 316.10 317.01
## 1962 316.69 317.69
## 1963 317.12 318.31
## 1964 317.79 318.71
## 1965 318.87 319.42
## 1966 319.79 321.08
## 1967 320.72 321.96
## 1968 321.31 322.84
## 1969 322.85 324.11
## 1970 323.98 325.13
## 1971 324.80 326.01
## 1972 326.50 327.55
## 1973 328.16 328.64
## 1974 328.30 329.58
## 1975 329.37 330.58
## 1976 330.15 331.62
## 1977 332.41 333.60
## 1978 333.75 334.90
## 1979 335.31 336.81
## 1980 337.27 338.32
## 1981 338.58 339.88
## 1982 339.60 340.90
## 1983 341.59 343.05
## 1984 343.10 344.70
## 1985 344.49 345.88
## 1986 345.93 347.22
## 1987 347.96 349.16
## 1988 350.20 351.41
## 1989 351.49 352.85
## 1990 353.07 354.33
## 1991 353.81 355.12
## 1992 354.32 355.57
## 1993 355.41 356.91
## 1994 357.53 358.98
## 1995 359.45 360.68
## 1996 360.89 362.24
## 1997 362.49 364.38
## 1998 365.52 367.13
## 1999 366.72 368.10
## 2000 368.37 369.67
## 2001 369.76 371.24
## 2002 372.25 373.79
## 2003 374.66 375.99
## 2004 376.16 377.51
## 2005 378.35 380.11
## 2006 380.15 381.82
## 2007 382.45 383.94
## 2008 384.19 385.56
## 2009 386.02 387.42
## 2010 388.67 389.79
## 2011 390.28 391.86
## 2012 392.98 394.34
## 2013 395.16 396.84
## 2014 397.28 398.91
## 2015 400.16 401.85
## 2016 403.55 404.45
## 2017 405.12 406.81
## 2018 408.02 409.07
## 2019 410.25 411.76
## 2020# Jan Feb Mar Apr May Jun Jul Aug
#1958 315.71 317.45 317.50 317.10 315.86 314.93
#1959 315.62 316.38 316.71 317.72 318.29 318.15 316.54 314.80
#1960 316.43 316.97 317.58 319.02 320.03 319.59 318.18 315.91
#Decompose a time series into components of trend,
#seasonal cycle, and random residual
co2.decompose <- decompose(co2.ts)
#Plot the time series and its three components
plot(co2.decompose, xlab ="")
#Change the title of the plot
my_plot.decomposed.ts = function(x, title="", ...) {
xx <- x$x
if (is.null(xx))
xx <- with(x, if (type == "additive")
random + trend + seasonal
else random * trend * seasonal)
plot(cbind(observed = xx, trend = x$trend,
seasonal = x$seasonal, random = x$random),
main=title, ...)
}
par(mar=c(2,4.5,2.5,0.2))
my_plot.decomposed.ts(co2.decompose,
expression(paste("Monthly Atmospheric ",
CO[2]," [ppm] at Mauna Loa Observatory")),
xlab ="", cex.lab=1.5, cex.axis=1.8)
plot(co2.decompose$trend, type = 'l')
lines(mon, co2m[,6], col = 'red')
plot(co2.decompose$trend - co2m[,6], type = 'l')
#Trend component
summary(co2.decompose$trend, na.rm=TRUE)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 315.4 330.1 352.9 355.2 377.9 412.8 12# Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
#315.4 330.1 352.9 355.2 377.9 412.8 12
plot(co2.decompose$trend, type = 'l')
#Seasonal component
round(summary(co2.decompose$seasonal, na.rm=TRUE),
digits = 3)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -3.256 -1.484 0.678 0.013 2.318 3.031#Min. 1st Qu. Median Mean 3rd Qu. Max.
#-3.256 -1.484 0.678 0.013 2.318 3.031
plot(co2.decompose$seasonal, type = 'l')
plot(co2.decompose$seasonal[1:50], type = 'l')
#Random component
round(summary(co2.decompose$random, na.rm=TRUE),
digits = 4)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -0.9628 -0.1817 0.0034 -0.0009 0.1818 1.2659 12# Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
#-0.9628 -0.1817 0.0034 -0.0009 0.1818 1.2659 12
which (co2.decompose$random > 1.0)## [1] 186 698#[1] 186 698
co2m[c(186, 698),]## year mon date average interpolated trend days
## 186 1973 8 1973.625 329.31 329.31 330.42 -1
## 698 2016 4 2016.292 407.45 407.45 404.60 25#year mon date average interpolated trend days
#186 1973 8 1973.625 329.31 329.31 330.42 -1
#698 2016 4 2016.292 407.45 407.45 404.60 25
sd(co2.decompose$random, na.rm=TRUE)## [1] 0.2953659#[1] 0.2953659
plot(co2.decompose$random, type = 'l')
###
min(co2.decompose$trend - co2m[,6], na.rm=TRUE)## [1] -0.9975###older version to be deleted
#setwd("/Users/sshen/climstats")
co2m = read.table("data/co2m.txt", header = TRUE)
dim(co2m)## [1] 749 7# [1] 749 7
co2m[1:3,]## year mon date average interpolated trend days
## 1 1958 3 1958.208 315.71 315.71 314.62 -1
## 2 1958 4 1958.292 317.45 317.45 315.29 -1
## 3 1958 5 1958.375 317.50 317.50 314.71 -1#year mon date average interpolated trend days
#1 1958 3 1958.208 315.71 315.71 314.62 -1
#2 1958 4 1958.292 317.45 317.45 315.29 -1
#3 1958 5 1958.375 317.50 317.50 314.71 -1
mon = co2m[,3]
co2 = co2m[,5]
##### test
co2_12 = c(1:2, co2, 8:12)
length(co2_12)/12## [1] 63#[1] 63
mon[35: (35 + 12*50 -1)] #Jan 1961-Dec 2010## [1] 1961.042 1961.125 1961.208 1961.292 1961.375 1961.458 1961.542 1961.625
## [9] 1961.708 1961.792 1961.875 1961.958 1962.042 1962.125 1962.208 1962.292
## [17] 1962.375 1962.458 1962.542 1962.625 1962.708 1962.792 1962.875 1962.958
## [25] 1963.042 1963.125 1963.208 1963.292 1963.375 1963.458 1963.542 1963.625
## [33] 1963.708 1963.792 1963.875 1963.958 1964.042 1964.125 1964.208 1964.292
## [41] 1964.375 1964.458 1964.542 1964.625 1964.708 1964.792 1964.875 1964.958
## [49] 1965.042 1965.125 1965.208 1965.292 1965.375 1965.458 1965.542 1965.625
## [57] 1965.708 1965.792 1965.875 1965.958 1966.042 1966.125 1966.208 1966.292
## [65] 1966.375 1966.458 1966.542 1966.625 1966.708 1966.792 1966.875 1966.958
## [73] 1967.042 1967.125 1967.208 1967.292 1967.375 1967.458 1967.542 1967.625
## [81] 1967.708 1967.792 1967.875 1967.958 1968.042 1968.125 1968.208 1968.292
## [89] 1968.375 1968.458 1968.542 1968.625 1968.708 1968.792 1968.875 1968.958
## [97] 1969.042 1969.125 1969.208 1969.292 1969.375 1969.458 1969.542 1969.625
## [105] 1969.708 1969.792 1969.875 1969.958 1970.042 1970.125 1970.208 1970.292
## [113] 1970.375 1970.458 1970.542 1970.625 1970.708 1970.792 1970.875 1970.958
## [121] 1971.042 1971.125 1971.208 1971.292 1971.375 1971.458 1971.542 1971.625
## [129] 1971.708 1971.792 1971.875 1971.958 1972.042 1972.125 1972.208 1972.292
## [137] 1972.375 1972.458 1972.542 1972.625 1972.708 1972.792 1972.875 1972.958
## [145] 1973.042 1973.125 1973.208 1973.292 1973.375 1973.458 1973.542 1973.625
## [153] 1973.708 1973.792 1973.875 1973.958 1974.042 1974.125 1974.208 1974.292
## [161] 1974.375 1974.458 1974.542 1974.625 1974.708 1974.792 1974.875 1974.958
## [169] 1975.042 1975.125 1975.208 1975.292 1975.375 1975.458 1975.542 1975.625
## [177] 1975.708 1975.792 1975.875 1975.958 1976.042 1976.125 1976.208 1976.292
## [185] 1976.375 1976.458 1976.542 1976.625 1976.708 1976.792 1976.875 1976.958
## [193] 1977.042 1977.125 1977.208 1977.292 1977.375 1977.458 1977.542 1977.625
## [201] 1977.708 1977.792 1977.875 1977.958 1978.042 1978.125 1978.208 1978.292
## [209] 1978.375 1978.458 1978.542 1978.625 1978.708 1978.792 1978.875 1978.958
## [217] 1979.042 1979.125 1979.208 1979.292 1979.375 1979.458 1979.542 1979.625
## [225] 1979.708 1979.792 1979.875 1979.958 1980.042 1980.125 1980.208 1980.292
## [233] 1980.375 1980.458 1980.542 1980.625 1980.708 1980.792 1980.875 1980.958
## [241] 1981.042 1981.125 1981.208 1981.292 1981.375 1981.458 1981.542 1981.625
## [249] 1981.708 1981.792 1981.875 1981.958 1982.042 1982.125 1982.208 1982.292
## [257] 1982.375 1982.458 1982.542 1982.625 1982.708 1982.792 1982.875 1982.958
## [265] 1983.042 1983.125 1983.208 1983.292 1983.375 1983.458 1983.542 1983.625
## [273] 1983.708 1983.792 1983.875 1983.958 1984.042 1984.125 1984.208 1984.292
## [281] 1984.375 1984.458 1984.542 1984.625 1984.708 1984.792 1984.875 1984.958
## [289] 1985.042 1985.125 1985.208 1985.292 1985.375 1985.458 1985.542 1985.625
## [297] 1985.708 1985.792 1985.875 1985.958 1986.042 1986.125 1986.208 1986.292
## [305] 1986.375 1986.458 1986.542 1986.625 1986.708 1986.792 1986.875 1986.958
## [313] 1987.042 1987.125 1987.208 1987.292 1987.375 1987.458 1987.542 1987.625
## [321] 1987.708 1987.792 1987.875 1987.958 1988.042 1988.125 1988.208 1988.292
## [329] 1988.375 1988.458 1988.542 1988.625 1988.708 1988.792 1988.875 1988.958
## [337] 1989.042 1989.125 1989.208 1989.292 1989.375 1989.458 1989.542 1989.625
## [345] 1989.708 1989.792 1989.875 1989.958 1990.042 1990.125 1990.208 1990.292
## [353] 1990.375 1990.458 1990.542 1990.625 1990.708 1990.792 1990.875 1990.958
## [361] 1991.042 1991.125 1991.208 1991.292 1991.375 1991.458 1991.542 1991.625
## [369] 1991.708 1991.792 1991.875 1991.958 1992.042 1992.125 1992.208 1992.292
## [377] 1992.375 1992.458 1992.542 1992.625 1992.708 1992.792 1992.875 1992.958
## [385] 1993.042 1993.125 1993.208 1993.292 1993.375 1993.458 1993.542 1993.625
## [393] 1993.708 1993.792 1993.875 1993.958 1994.042 1994.125 1994.208 1994.292
## [401] 1994.375 1994.458 1994.542 1994.625 1994.708 1994.792 1994.875 1994.958
## [409] 1995.042 1995.125 1995.208 1995.292 1995.375 1995.458 1995.542 1995.625
## [417] 1995.708 1995.792 1995.875 1995.958 1996.042 1996.125 1996.208 1996.292
## [425] 1996.375 1996.458 1996.542 1996.625 1996.708 1996.792 1996.875 1996.958
## [433] 1997.042 1997.125 1997.208 1997.292 1997.375 1997.458 1997.542 1997.625
## [441] 1997.708 1997.792 1997.875 1997.958 1998.042 1998.125 1998.208 1998.292
## [449] 1998.375 1998.458 1998.542 1998.625 1998.708 1998.792 1998.875 1998.958
## [457] 1999.042 1999.125 1999.208 1999.292 1999.375 1999.458 1999.542 1999.625
## [465] 1999.708 1999.792 1999.875 1999.958 2000.042 2000.125 2000.208 2000.292
## [473] 2000.375 2000.458 2000.542 2000.625 2000.708 2000.792 2000.875 2000.958
## [481] 2001.042 2001.125 2001.208 2001.292 2001.375 2001.458 2001.542 2001.625
## [489] 2001.708 2001.792 2001.875 2001.958 2002.042 2002.125 2002.208 2002.292
## [497] 2002.375 2002.458 2002.542 2002.625 2002.708 2002.792 2002.875 2002.958
## [505] 2003.042 2003.125 2003.208 2003.292 2003.375 2003.458 2003.542 2003.625
## [513] 2003.708 2003.792 2003.875 2003.958 2004.042 2004.125 2004.208 2004.292
## [521] 2004.375 2004.458 2004.542 2004.625 2004.708 2004.792 2004.875 2004.958
## [529] 2005.042 2005.125 2005.208 2005.292 2005.375 2005.458 2005.542 2005.625
## [537] 2005.708 2005.792 2005.875 2005.958 2006.042 2006.125 2006.208 2006.292
## [545] 2006.375 2006.458 2006.542 2006.625 2006.708 2006.792 2006.875 2006.958
## [553] 2007.042 2007.125 2007.208 2007.292 2007.375 2007.458 2007.542 2007.625
## [561] 2007.708 2007.792 2007.875 2007.958 2008.042 2008.125 2008.208 2008.292
## [569] 2008.375 2008.458 2008.542 2008.625 2008.708 2008.792 2008.875 2008.958
## [577] 2009.042 2009.125 2009.208 2009.292 2009.375 2009.458 2009.542 2009.625
## [585] 2009.708 2009.792 2009.875 2009.958 2010.042 2010.125 2010.208 2010.292
## [593] 2010.375 2010.458 2010.542 2010.625 2010.708 2010.792 2010.875 2010.958co2mat = matrix(co2_12, ncol =12, byrow = TRUE)
clim = matrix(0, nrow =length(mon), ncol = 12)
clim = colMeans(co2mat)
plot(clim)
co2trend = 1:length(mon)
for (i in 1: length(mon)){
co2trend[i] = co2[i]-clim[co2m[i,2]] + mean(clim)
}
plot(co2trend, type = 'l')
lines(co2m[,6], type = 'l', col ='red')
plot(co2trend - co2m[,6], type = 'l')
plot(co2m[,6], type = 'l')
####
#Plot the nonlinear trend line
#install.packages('forecast')
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
co2.ts = ts(co2, start=c(1958,3), end=c(2020,7),
frequency =12)
co2.decompose <- decompose(co2.ts)
plot(co2.decompose$trend,
ylab = 'CO2 [ppm]', xlab = 'Time [mon]')
R Plot Fig. 7.3: Time Series Forecast for CO2
#by ETS(Error, Trend, and Seasonal) smoothing
library(forecast)
co2forecasts <- ets(co2.ts,model="AAA")
#ETS model AAA means exponential smoothing with
#season and trend
#Forecast 48 months to the future
co2forecasts <- forecast(co2forecasts, h=48)
plot(co2forecasts,
ylab="Parts per million [ppm]",
main = "Observed (red) and ETS Forecast (blue)",
cex.axis = 1.5, cex.lab = 1.5, cex.main = 1.2,
lwd = 2)
lines(co2.ts, col='red')
R Plot Fig. 7.4: Daily Tmin at ST PAUL Station, Minnesota
#1941-1950 MINNEAPOLIS/ST PAUL
#GHCND:USW00014927 44.8831N 93.2289W 265.8Elev
#setwd("/Users/sshen/climstats")
tda = read.table("data/StPaulStn.txt", header = F)
dim(tda)## [1] 7589 13#[1] 7589 13
#Col12 = Tmin, col9 = PRCP, col8 = Date
t1 = tda[,8] #date
tmin = tda[, 12] #Tmin data
t1[1462]## [1] 19410101#[1] 19410101
t1[5112]## [1] 19501231#[1] 19501231
da = tmin[1462:5112]
ta = t1[1462:5112]
dav = as.vector(t(da))
Tmin.ts = ts(dav, start=c(1941,1,1),
end=c(1950,12,31),
frequency =365)
#ETS time series decomposition
Tmin.decompose <- decompose(Tmin.ts)
plot(Tmin.decompose)
#Change the title of the ETS decomposition figure
my_plot.decomposed.ts = function(x, title="", ...) {
xx <- x$x
if (is.null(xx))
xx <- with(x, if (type == "additive")
random + trend + seasonal
else random * trend * seasonal)
plot(cbind(observed = xx, trend = x$trend,
seasonal = x$seasonal, random = x$random),
main=title, ...)
}
par(mar=c(2,4.5,2.5,0.2))
my_plot.decomposed.ts(Tmin.decompose,
title = "St. Paul (Minnesota, USA) Daily Tmin (deg C): 1941-1950",
xlab ="", cex.lab=1.5, cex.axis=1.8)
#Check the minimum noise day:
n1 = which(Tmin.decompose$random < -20)
tda[(n1 + 1641 - 3): (n1 + 1641 + 3), 8:13]## V8 V9 V10 V11 V12 V13
## 4265 19480904 0.0 -9999 32.8 17.2 -9999
## 4266 19480905 0.0 -9999 32.2 17.2 -9999
## 4267 19480906 0.0 -9999 25.6 16.7 -9999
## 4268 19480907 0.0 -9999 23.3 13.3 -9999
## 4269 19480908 3.6 -9999 24.4 12.8 -9999
## 4270 19480909 0.0 -9999 21.1 10.0 -9999
## 4271 19480910 0.0 -9999 22.2 6.1 -9999#4268 19480907 0 -9999 23.3 13.3 -9999
#PRCP = 0 [mm], Tmin = 13.3 [deg C]
R Plot Fig. 7.5: Generate White Noise
set.seed(119) # set seed to ensure the same simulation result
wa <- rnorm(500, mean = 0, sd = 1)
par(mar=c(4.5,4.5,2.5, 0.3))
plot(wa, type = 'l',
xlab = 'Time', ylab = expression(W[t]),
main = 'A Simulated Time Series of White Noise ',
cex.lab=1.5, cex.axis=1.5)
R Plot Fig. 7.6: Histogram and ACF of White Noise
par(mfrow=c(1,2))
par(mar=c(4.5,4.5,2.5, 0))
hist(wa, freq=F, breaks=16,
xlim =c(-4,4), ylim=c(0,0.4),
xlab = expression(W[t]),
main='Histogram of a White Noise Time Series',
cex.lab=1.5, cex.axis=1.5)
x=seq(-4,4, by=0.1)
lines(x,dnorm(x, mean=0, sd=1), col='blue',
type = 'l', lwd=1.5)
par(mar=c(4.5,4.5,3, 0.3))
acf(wa, main='Auto-correlation of White Noise',
cex.lab=1.5, cex.axis=1.5)
#dev.off()
R Plot Fig. 7.7: Random Walk Time Series
#(a) Five realizations with different sigma values
n = 1000 #Total number of time steps
m = 5 #Number of time series realizations
a = 0.05 #Drift delta = 0.05
b = c(0.1, 0.4, 0.7, 1.0, 1.5) #SD sigma
#Generate the random walk time series data
x = matrix(rep(0, m*n), nrow=m)
for(i in 1:m){
w = rnorm(n, mean =0, sd = b[i])
for(j in 1:(n-1)){
x[i,j+1] = a + x[i,j] + w[j]
}
}
#Plot the five time series realizations
par(mfrow=c(1,2))
par(mar=c(4.5,4.5, 2.5, 0.5))
plot(x[1,], type='l', ylim=c(-20,100),
xlab ="Time steps: t", ylab = expression(X[t]),
main = expression('Five realizations of random walks:'
~ delta ~'= 0.05'),
cex.lab=1.6, cex.axis=1.5, cex.main =1.6)
lines(x[2,], type='l', col='blue')
lines(x[3,], type='l', col='red')
lines(x[4,], type='l', col='orange')
lines(x[5,], type='l', col='brown')
legend(-100, 110,
legend=c(expression(sigma ~ '= 0.1'),
expression(sigma ~ '= 0.4'),
expression(sigma ~ '= 0.7'),
expression(sigma ~ '= 1.0'),
expression(sigma ~ '= 1.5')),
col=c('black','blue','red','orange','brown'),
x.intersp = 0.2, y.intersp = 0.4,
seg.len = 0.6, lty =1, bty ='n', cex = 1.3)
text(20,-20, "(a)", cex =1.4)
#(b) 100 realizations with fixed parameters
#Random walk to show variance increasing with time
n = 1000 #Total number of time steps
m = 100 #Number of time series realizations
a = 0.0 #Drift delta = 0
b = rep(1, m) #SD sigma is the same
#Generate the random walk time series data
x = matrix(rep(0, m*n), nrow=m)
for(i in 1:m){
w = rnorm(n, mean =0, sd = b[i])
for(j in 1:(n-1)){
x[i,j+1] = a + x[i,j] + w[j]
}
}
#Plot the series realizations
par(mar=c(4.5,4.5, 2.5, 0.8))
plot(x[1,], type='l', ylim=c(-100,100),
xlab ="Time steps: t", ylab = expression(X[t]),
main = expression('100 realizations of random walks:'
~ delta ~'= 0,'~~ sigma ~'= 1'),
cex.lab=1.6, cex.axis=1.5, cex.main =1.6)
for(i in 2:m){
lines(x[i,], type='l')
}
library(matrixStats)
y = colSds(x)
lines(y,type='l', lwd=2, col='red')
lines(-y,type='l', lwd=2, col='red')
z = sqrt(1:n)
lines(z,type='l', lwd=2, col='blue')
lines(-z,type='l', lwd=2, col='blue')
legend(-150, 120,
legend=c('Standard deviation of the simulations SD(t)',
expression('Theorectical formula: SD(t)'%prop% sqrt(t))),
col=c('red','blue'), cex = 1.3,
x.intersp = 0.2, y.intersp = 0.3,
seg.len = 0.4, lty=1, bty='n', lwd =2)
text(20,-100, "(b)", cex =1.4)
#dev.off() #go back to R's default figure setting
R Plot Fig. 7.8: Moving Average Time Series
n = 121
w = rnorm(n)
x1 = x2 = x3 = rep(0, n)
for (t in 5:n){
x1[t] = 0.5*w[t] + 0.5*w[t-1]
x2[t] = 0.9*w[t] + 0.1*w[t-1]
x3[t] = (1/5)*(w[t] + w[t-1] + w[t-2] + w[t-3] + w[t-4])
}
plot.ts(x1[5:n], ylim= c(-2,3),
main="Realizations of moving average time series",
xlab ="Time steps", ylab = expression(X[t]),
cex.lab=1.5, cex.axis=1.5, lwd=2)
lines(x2[5:n], col='blue')
lines(x3[5:n], col='red', lwd=2)
legend(5, 3.2, bty='n', lty=1, y.intersp=0.7,
legend=(c("MA(1): a = 0.5, b = 0.5",
"MA(1): a = 0.9, b = 0.1",
"MA(5): Uniform weight = 1/5")),
col=c('black','blue','red'),
lwd = c(2,1,2))
R Plot Fig. 7.9: Autoregressive Time Series AR(1)
set.seed(791)
n = 121 #Number of time steps
m = 4 #Number of realizations
lam = 0.8 #Decay parameter
x = matrix(rep(4, n*m), nrow=m) #x0 = 4
#Simulate the time series data
for(k in 1:m){
for(i in 1:(n-1)){
x[k,i+1] = x[k,i]*lam +
rnorm(1, mean=0, sd=0.25)
}
}
#Plot the realizations and their mean
plot.ts(x[1,], type = 'l', ylim=c(-2,4),
main="Realizations of an AR(1) time series and their mean",
xlab ="Time steps", ylab = expression(X[t]),
cex.lab=1.5, cex.axis=1.5,
col='purple')
lines(x[2,], col='blue')
lines(x[3,], col='red')
lines(x[4,], col='orange')
lines(colMeans(x), lwd=3)
R Plot Fig. 7.10: ACF of AR(1)
#setwd('/Users/sshen/climstats')
n = 481 #Number of time steps
m = 2 #Number of realizations
lam = c(0.9, 0.6) #Decay parameter
x = matrix(rep(4, n*m), nrow=m) #x0 = 4
#Simulate the time series data
for(k in 1:m){
for(i in 1:(n-1)){
x[k,i+1] = x[k,i]*lam[k] +
rnorm(1, mean=0, sd=0.25)
}
}
#Plot the auto-correlation function
#setEPS() #Automatically saves the .eps file
#postscript("fig0710.eps", width=10, height=5)
par(mfrow=c(1,2))
par(mar=c(4.5, 4.5, 3, 1))
acf(x[1,], lag.max = 36,
main = 'Auto-correlation function of AR(1)',
xlab='Time lag',
cex.lab=1.5, cex.axis=1.5)
text(20, 0.8, bquote('Decay parameter'~lambda == 0.9),
cex=1.5)
par(mar=c(4.5,4.5,3, 0.3))
acf(x[2,], lag.max = 36,
main = 'Auto-correlation function of AR(1)',
xlab='Time lag', col='red',
cex.lab=1.5, cex.axis=1.5)
text(20, 0.8, expression('Decay parameter'~lambda == 0.6),
col='red', cex=1.5)
#dev.off()
R Plot Fig. 7.11: ARIMA Model Fitting
#setwd("/Users/sshen/climstats")
co2m = read.table("data/co2m.txt", header = TRUE)
mon = co2m[,3]
co2 = co2m[,5]
co2.ts = ts(co2, start=c(1958,3), end=c(2020,7),
frequency =12)
#Then fit an AR(1) model
co2.AR1 <- arima(co2.ts, order = c(1,0,0))
#Obtain the AR(1) model fit data
AR1_fit <- co2.ts - residuals(co2.AR1)
#Fit an MA(1) model
co2.MA1 <- arima(co2.ts, order = c(0,0,1))
#Obtain the MA(1) model fit data
MA1_fit <- co2.ts - residuals(co2.MA1)
#setEPS() #Automatically saves the .eps file
#postscript("fig0711.eps", width=8, height=6)
#Plot the CO2 time series
par(mar=c(3,4.5,2.5, 0.5))
plot(mon, co2, ylim=c(310,420),
type='l', lwd=2, col = 'red',
main = bquote("Monthly Atmospheric "~
CO[2]~" at Mauna Loa Observatory"),
xlab ="",
ylab="Parts Per Million [ppm]",
cex.main=1.5,
cex.lab =1.5, cex.axis=1.5)
#Plot the AR(1) model data on the CO2 data
points(AR1_fit, type = "l", lwd=1,
col = 'black', lty = 2)
points(MA1_fit, type = "l", lwd=1,
col = 'blue', lty = 2)
legend(1960, 420, lty=c(1, 2, 2),
col=c('red','black', 'blue'),
legend=c("Observed data",
"AR(1) model fitting",
"MA(1) model fitting"),
text.col = c('red','black', 'blue'),
x.intersp = 0.2, y.intersp = 1.4,
seg.len = 0.8, cex=1.5, bty = 'n')
#dev.off()