PM10 Dataset Example

This vignette shows how to use the built-in pm10 dataset with tsqn.

library(tsqn)
#> Loading required package: robustbase
#> Loading required package: MASS
#> Loading required package: fracdiff
data("pm10")

dim(pm10)
#> [1] 1826    8
head(pm10)
#>   Laranjeiras Carapina Camburi     Sua VixCentro    Ibes VVCentro Cariacica
#> 1     24.0000  14.5833 16.1250 20.6667   18.0833 15.5417  21.3333   31.1667
#> 2     21.7917  14.5000 22.5000 26.5417   21.5417 17.4167  17.5417   32.6250
#> 3     31.7083  19.3333 28.2917 29.0417   28.2917 28.9167  39.8333   50.5417
#> 4     24.5833  22.9583 21.3750 20.8750   23.7917 19.2917  26.9583   38.1250
#> 5     34.5417  19.5000 28.9583 31.9583   31.1667 21.3333  37.9167   46.7083
#> 6     37.0000  17.1667 26.1250 29.0833   33.3333 23.4167  38.2917   40.4167

The complete dataset has 1826 observations for 8 monitoring stations. For faster examples in this vignette, we use the first 365 observations.

pm10_subset <- as.matrix(pm10[1:365, ])

qn_cor <- corMatQn(pm10_subset)
qn_cov <- covMatQn(pm10_subset)

round(qn_cor, 3)
#>       [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]
#> [1,] 1.000 0.571 0.584 0.604 0.499 0.687 0.399 0.557
#> [2,] 0.571 1.000 0.720 0.654 0.667 0.705 0.431 0.757
#> [3,] 0.584 0.720 1.000 0.643 0.654 0.655 0.431 0.631
#> [4,] 0.604 0.654 0.643 1.000 0.759 0.689 0.476 0.630
#> [5,] 0.499 0.667 0.654 0.759 1.000 0.547 0.583 0.714
#> [6,] 0.687 0.705 0.655 0.689 0.547 1.000 0.301 0.709
#> [7,] 0.399 0.431 0.431 0.476 0.583 0.301 1.000 0.375
#> [8,] 0.557 0.757 0.631 0.630 0.714 0.709 0.375 1.000
round(qn_cov, 1)
#>       [,1] [,2] [,3] [,4] [,5] [,6]  [,7]  [,8]
#> [1,] 114.0 39.7 49.5 49.8 39.0 66.7  55.2  83.1
#> [2,]  39.7 39.2 33.8 31.9 28.9 40.3  36.8  64.1
#> [3,]  49.5 33.8 57.0 38.0 33.5 44.7  42.9  62.3
#> [4,]  49.8 31.9 38.0 55.6 40.2 46.8  48.8  65.7
#> [5,]  39.0 28.9 33.5 40.2 49.0 34.6  52.6  66.3
#> [6,]  66.7 40.3 44.7 46.8 34.6 78.1  35.0  85.9
#> [7,]  55.2 36.8 42.9 48.8 52.6 35.0 136.7  69.2
#> [8,]  83.1 64.1 62.3 65.7 66.3 85.9  69.2 168.5

Robust ACF and robust spectral analysis for one station:

vix <- pm10_subset[, "VixCentro"]

acf_qn <- robacf(vix, lag.max = 24, type = "correlation", plot = FALSE)
head(acf_qn$acf[, 1, 1], 10)
#>  [1]  1.000000000  0.272263665  0.009576851  0.095893596  0.055512722
#>  [6] -0.102051424 -0.075357148  0.036692338 -0.037377451 -0.009397082

per_qn <- PerQn(vix)
length(per_qn)
#> [1] 363
head(per_qn, 10)
#>  [1] 10.8263433  9.2346587  9.9512860 14.0949867 19.2784413 21.1307737
#>  [7] 16.9945525  8.6261808  1.4026843  0.4378287

GPH_estimate(vix, method = "GPH-Qn")
#> $method
#> [1] "Qn"
#> 
#> $d
#> [1] 0.04464066
#> 
#> $sd.reg
#> [1] 0.06856734
#> 
#> [[4]]
#> [1] 0.7
tsqn:::plot.robacf(acf_qn, main = "PM10 Robust ACF (VixCentro)")

For comparison, the classical (standard) ACF is:

stats::acf(vix, lag.max = 24, main = "PM10 Standard ACF (VixCentro)")