R version 2.13.0 (2011-04-13)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(16111,15554,15220,14807,14291,14653,17006,18032,16558,16102,15055,15484,14596,14609,13923,14226,14056,14278,16142,16509,15680,14086,13129,13086,13096,12280,11534,11135,10903,10926,13220,13581,11788,11088,10434,11061,10828,10270,10360,9899,9395,9944,12117,12474,11106,10643,10227,11273,11516,11583,11605,11414,11181,12000,14007,14582,13251,12806,12645,13869,13342,13079,12513,12331,11882,12388,14394,14635,13218,12554,12031,12429),dim=c(1,72),dimnames=list(c('Werkloosheid'),1:72))
> y <- array(NA,dim=c(1,72),dimnames=list(c('Werkloosheid'),1:72))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> library(lattice)
> library(lmtest)
Loading required package: zoo
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Werkloosheid M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 16111 1 0 0 0 0 0 0 0 0 0 0
2 15554 0 1 0 0 0 0 0 0 0 0 0
3 15220 0 0 1 0 0 0 0 0 0 0 0
4 14807 0 0 0 1 0 0 0 0 0 0 0
5 14291 0 0 0 0 1 0 0 0 0 0 0
6 14653 0 0 0 0 0 1 0 0 0 0 0
7 17006 0 0 0 0 0 0 1 0 0 0 0
8 18032 0 0 0 0 0 0 0 1 0 0 0
9 16558 0 0 0 0 0 0 0 0 1 0 0
10 16102 0 0 0 0 0 0 0 0 0 1 0
11 15055 0 0 0 0 0 0 0 0 0 0 1
12 15484 0 0 0 0 0 0 0 0 0 0 0
13 14596 1 0 0 0 0 0 0 0 0 0 0
14 14609 0 1 0 0 0 0 0 0 0 0 0
15 13923 0 0 1 0 0 0 0 0 0 0 0
16 14226 0 0 0 1 0 0 0 0 0 0 0
17 14056 0 0 0 0 1 0 0 0 0 0 0
18 14278 0 0 0 0 0 1 0 0 0 0 0
19 16142 0 0 0 0 0 0 1 0 0 0 0
20 16509 0 0 0 0 0 0 0 1 0 0 0
21 15680 0 0 0 0 0 0 0 0 1 0 0
22 14086 0 0 0 0 0 0 0 0 0 1 0
23 13129 0 0 0 0 0 0 0 0 0 0 1
24 13086 0 0 0 0 0 0 0 0 0 0 0
25 13096 1 0 0 0 0 0 0 0 0 0 0
26 12280 0 1 0 0 0 0 0 0 0 0 0
27 11534 0 0 1 0 0 0 0 0 0 0 0
28 11135 0 0 0 1 0 0 0 0 0 0 0
29 10903 0 0 0 0 1 0 0 0 0 0 0
30 10926 0 0 0 0 0 1 0 0 0 0 0
31 13220 0 0 0 0 0 0 1 0 0 0 0
32 13581 0 0 0 0 0 0 0 1 0 0 0
33 11788 0 0 0 0 0 0 0 0 1 0 0
34 11088 0 0 0 0 0 0 0 0 0 1 0
35 10434 0 0 0 0 0 0 0 0 0 0 1
36 11061 0 0 0 0 0 0 0 0 0 0 0
37 10828 1 0 0 0 0 0 0 0 0 0 0
38 10270 0 1 0 0 0 0 0 0 0 0 0
39 10360 0 0 1 0 0 0 0 0 0 0 0
40 9899 0 0 0 1 0 0 0 0 0 0 0
41 9395 0 0 0 0 1 0 0 0 0 0 0
42 9944 0 0 0 0 0 1 0 0 0 0 0
43 12117 0 0 0 0 0 0 1 0 0 0 0
44 12474 0 0 0 0 0 0 0 1 0 0 0
45 11106 0 0 0 0 0 0 0 0 1 0 0
46 10643 0 0 0 0 0 0 0 0 0 1 0
47 10227 0 0 0 0 0 0 0 0 0 0 1
48 11273 0 0 0 0 0 0 0 0 0 0 0
49 11516 1 0 0 0 0 0 0 0 0 0 0
50 11583 0 1 0 0 0 0 0 0 0 0 0
51 11605 0 0 1 0 0 0 0 0 0 0 0
52 11414 0 0 0 1 0 0 0 0 0 0 0
53 11181 0 0 0 0 1 0 0 0 0 0 0
54 12000 0 0 0 0 0 1 0 0 0 0 0
55 14007 0 0 0 0 0 0 1 0 0 0 0
56 14582 0 0 0 0 0 0 0 1 0 0 0
57 13251 0 0 0 0 0 0 0 0 1 0 0
58 12806 0 0 0 0 0 0 0 0 0 1 0
59 12645 0 0 0 0 0 0 0 0 0 0 1
60 13869 0 0 0 0 0 0 0 0 0 0 0
61 13342 1 0 0 0 0 0 0 0 0 0 0
62 13079 0 1 0 0 0 0 0 0 0 0 0
63 12513 0 0 1 0 0 0 0 0 0 0 0
64 12331 0 0 0 1 0 0 0 0 0 0 0
65 11882 0 0 0 0 1 0 0 0 0 0 0
66 12388 0 0 0 0 0 1 0 0 0 0 0
67 14394 0 0 0 0 0 0 1 0 0 0 0
68 14635 0 0 0 0 0 0 0 1 0 0 0
69 13218 0 0 0 0 0 0 0 0 1 0 0
70 12554 0 0 0 0 0 0 0 0 0 1 0
71 12031 0 0 0 0 0 0 0 0 0 0 1
72 12429 0 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
12867.00 381.17 28.83 -341.17 -565.00 -915.67
M6 M7 M8 M9 M10 M11
-502.17 1614.00 2101.83 733.17 12.83 -613.50
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2625.8 -1400.6 -274.2 1570.4 3222.2
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12867.00 775.52 16.591 <2e-16 ***
M1 381.17 1096.76 0.348 0.7294
M2 28.83 1096.76 0.026 0.9791
M3 -341.17 1096.76 -0.311 0.7568
M4 -565.00 1096.76 -0.515 0.6083
M5 -915.67 1096.76 -0.835 0.4071
M6 -502.17 1096.76 -0.458 0.6487
M7 1614.00 1096.76 1.472 0.1464
M8 2101.83 1096.76 1.916 0.0601 .
M9 733.17 1096.76 0.668 0.5064
M10 12.83 1096.76 0.012 0.9907
M11 -613.50 1096.76 -0.559 0.5780
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1900 on 60 degrees of freedom
Multiple R-squared: 0.2049, Adjusted R-squared: 0.05909
F-statistic: 1.405 on 11 and 60 DF, p-value: 0.1943
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 0.23821260 4.764252e-01 7.617874e-01
[2,] 0.14404618 2.880924e-01 8.559538e-01
[3,] 0.08333404 1.666681e-01 9.166660e-01
[4,] 0.04988127 9.976255e-02 9.501187e-01
[5,] 0.03979482 7.958963e-02 9.602052e-01
[6,] 0.06202432 1.240486e-01 9.379757e-01
[7,] 0.06950026 1.390005e-01 9.304997e-01
[8,] 0.15206198 3.041240e-01 8.479380e-01
[9,] 0.22688192 4.537638e-01 7.731181e-01
[10,] 0.34243134 6.848627e-01 6.575687e-01
[11,] 0.47883247 9.576649e-01 5.211675e-01
[12,] 0.65361508 6.927698e-01 3.463849e-01
[13,] 0.78348282 4.330344e-01 2.165172e-01
[14,] 0.87851940 2.429612e-01 1.214806e-01
[15,] 0.92097349 1.580530e-01 7.902651e-02
[16,] 0.95147719 9.704563e-02 4.852281e-02
[17,] 0.96367838 7.264323e-02 3.632162e-02
[18,] 0.97426560 5.146879e-02 2.573440e-02
[19,] 0.98529805 2.940391e-02 1.470195e-02
[20,] 0.98948655 2.102690e-02 1.051345e-02
[21,] 0.99141935 1.716129e-02 8.580646e-03
[22,] 0.99233958 1.532085e-02 7.660423e-03
[23,] 0.99437881 1.124238e-02 5.621190e-03
[24,] 0.99658148 6.837037e-03 3.418518e-03
[25,] 0.99692892 6.142165e-03 3.071082e-03
[26,] 0.99771654 4.566923e-03 2.283462e-03
[27,] 0.99853063 2.938740e-03 1.469370e-03
[28,] 0.99908869 1.822616e-03 9.113082e-04
[29,] 0.99938078 1.238430e-03 6.192152e-04
[30,] 0.99963722 7.255502e-04 3.627751e-04
[31,] 0.99980678 3.864494e-04 1.932247e-04
[32,] 0.99989327 2.134697e-04 1.067349e-04
[33,] 0.99995832 8.336524e-05 4.168262e-05
[34,] 0.99998263 3.473333e-05 1.736667e-05
[35,] 0.99999142 1.715780e-05 8.578899e-06
[36,] 0.99999390 1.219978e-05 6.099888e-06
[37,] 0.99998658 2.684898e-05 1.342449e-05
[38,] 0.99997207 5.585117e-05 2.792559e-05
[39,] 0.99991511 1.697718e-04 8.488591e-05
[40,] 0.99961491 7.701701e-04 3.850851e-04
[41,] 0.99831464 3.370713e-03 1.685357e-03
[42,] 0.99206050 1.587900e-02 7.939499e-03
[43,] 0.96632048 6.735903e-02 3.367952e-02
> postscript(file="/var/wessaorg/rcomp/tmp/1e32d1322553197.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/2tkwi1322553197.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/3f0h21322553197.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/4o48y1322553197.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/5nnai1322553197.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 72
Frequency = 1
1 2 3 4 5 6
2862.83333 2658.16667 2694.16667 2505.00000 2339.66667 2288.16667
7 8 9 10 11 12
2525.00000 3063.16667 2957.83333 3222.16667 2801.50000 2617.00000
13 14 15 16 17 18
1347.83333 1713.16667 1397.16667 1924.00000 2104.66667 1913.16667
19 20 21 22 23 24
1661.00000 1540.16667 2079.83333 1206.16667 875.50000 219.00000
25 26 27 28 29 30
-152.16667 -615.83333 -991.83333 -1167.00000 -1048.33333 -1438.83333
31 32 33 34 35 36
-1261.00000 -1387.83333 -1812.16667 -1791.83333 -1819.50000 -1806.00000
37 38 39 40 41 42
-2420.16667 -2625.83333 -2165.83333 -2403.00000 -2556.33333 -2420.83333
43 44 45 46 47 48
-2364.00000 -2494.83333 -2494.16667 -2236.83333 -2026.50000 -1594.00000
49 50 51 52 53 54
-1732.16667 -1312.83333 -920.83333 -888.00000 -770.33333 -364.83333
55 56 57 58 59 60
-474.00000 -386.83333 -349.16667 -73.83333 391.50000 1002.00000
61 62 63 64 65 66
93.83333 183.16667 -12.83333 29.00000 -69.33333 23.16667
67 68 69 70 71 72
-87.00000 -333.83333 -382.16667 -325.83333 -222.50000 -438.00000
> postscript(file="/var/wessaorg/rcomp/tmp/6jzzq1322553197.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 72
Frequency = 1
lag(myerror, k = 1) myerror
0 2862.83333 NA
1 2658.16667 2862.83333
2 2694.16667 2658.16667
3 2505.00000 2694.16667
4 2339.66667 2505.00000
5 2288.16667 2339.66667
6 2525.00000 2288.16667
7 3063.16667 2525.00000
8 2957.83333 3063.16667
9 3222.16667 2957.83333
10 2801.50000 3222.16667
11 2617.00000 2801.50000
12 1347.83333 2617.00000
13 1713.16667 1347.83333
14 1397.16667 1713.16667
15 1924.00000 1397.16667
16 2104.66667 1924.00000
17 1913.16667 2104.66667
18 1661.00000 1913.16667
19 1540.16667 1661.00000
20 2079.83333 1540.16667
21 1206.16667 2079.83333
22 875.50000 1206.16667
23 219.00000 875.50000
24 -152.16667 219.00000
25 -615.83333 -152.16667
26 -991.83333 -615.83333
27 -1167.00000 -991.83333
28 -1048.33333 -1167.00000
29 -1438.83333 -1048.33333
30 -1261.00000 -1438.83333
31 -1387.83333 -1261.00000
32 -1812.16667 -1387.83333
33 -1791.83333 -1812.16667
34 -1819.50000 -1791.83333
35 -1806.00000 -1819.50000
36 -2420.16667 -1806.00000
37 -2625.83333 -2420.16667
38 -2165.83333 -2625.83333
39 -2403.00000 -2165.83333
40 -2556.33333 -2403.00000
41 -2420.83333 -2556.33333
42 -2364.00000 -2420.83333
43 -2494.83333 -2364.00000
44 -2494.16667 -2494.83333
45 -2236.83333 -2494.16667
46 -2026.50000 -2236.83333
47 -1594.00000 -2026.50000
48 -1732.16667 -1594.00000
49 -1312.83333 -1732.16667
50 -920.83333 -1312.83333
51 -888.00000 -920.83333
52 -770.33333 -888.00000
53 -364.83333 -770.33333
54 -474.00000 -364.83333
55 -386.83333 -474.00000
56 -349.16667 -386.83333
57 -73.83333 -349.16667
58 391.50000 -73.83333
59 1002.00000 391.50000
60 93.83333 1002.00000
61 183.16667 93.83333
62 -12.83333 183.16667
63 29.00000 -12.83333
64 -69.33333 29.00000
65 23.16667 -69.33333
66 -87.00000 23.16667
67 -333.83333 -87.00000
68 -382.16667 -333.83333
69 -325.83333 -382.16667
70 -222.50000 -325.83333
71 -438.00000 -222.50000
72 NA -438.00000
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 2658.16667 2862.83333
[2,] 2694.16667 2658.16667
[3,] 2505.00000 2694.16667
[4,] 2339.66667 2505.00000
[5,] 2288.16667 2339.66667
[6,] 2525.00000 2288.16667
[7,] 3063.16667 2525.00000
[8,] 2957.83333 3063.16667
[9,] 3222.16667 2957.83333
[10,] 2801.50000 3222.16667
[11,] 2617.00000 2801.50000
[12,] 1347.83333 2617.00000
[13,] 1713.16667 1347.83333
[14,] 1397.16667 1713.16667
[15,] 1924.00000 1397.16667
[16,] 2104.66667 1924.00000
[17,] 1913.16667 2104.66667
[18,] 1661.00000 1913.16667
[19,] 1540.16667 1661.00000
[20,] 2079.83333 1540.16667
[21,] 1206.16667 2079.83333
[22,] 875.50000 1206.16667
[23,] 219.00000 875.50000
[24,] -152.16667 219.00000
[25,] -615.83333 -152.16667
[26,] -991.83333 -615.83333
[27,] -1167.00000 -991.83333
[28,] -1048.33333 -1167.00000
[29,] -1438.83333 -1048.33333
[30,] -1261.00000 -1438.83333
[31,] -1387.83333 -1261.00000
[32,] -1812.16667 -1387.83333
[33,] -1791.83333 -1812.16667
[34,] -1819.50000 -1791.83333
[35,] -1806.00000 -1819.50000
[36,] -2420.16667 -1806.00000
[37,] -2625.83333 -2420.16667
[38,] -2165.83333 -2625.83333
[39,] -2403.00000 -2165.83333
[40,] -2556.33333 -2403.00000
[41,] -2420.83333 -2556.33333
[42,] -2364.00000 -2420.83333
[43,] -2494.83333 -2364.00000
[44,] -2494.16667 -2494.83333
[45,] -2236.83333 -2494.16667
[46,] -2026.50000 -2236.83333
[47,] -1594.00000 -2026.50000
[48,] -1732.16667 -1594.00000
[49,] -1312.83333 -1732.16667
[50,] -920.83333 -1312.83333
[51,] -888.00000 -920.83333
[52,] -770.33333 -888.00000
[53,] -364.83333 -770.33333
[54,] -474.00000 -364.83333
[55,] -386.83333 -474.00000
[56,] -349.16667 -386.83333
[57,] -73.83333 -349.16667
[58,] 391.50000 -73.83333
[59,] 1002.00000 391.50000
[60,] 93.83333 1002.00000
[61,] 183.16667 93.83333
[62,] -12.83333 183.16667
[63,] 29.00000 -12.83333
[64,] -69.33333 29.00000
[65,] 23.16667 -69.33333
[66,] -87.00000 23.16667
[67,] -333.83333 -87.00000
[68,] -382.16667 -333.83333
[69,] -325.83333 -382.16667
[70,] -222.50000 -325.83333
[71,] -438.00000 -222.50000
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 2658.16667 2862.83333
2 2694.16667 2658.16667
3 2505.00000 2694.16667
4 2339.66667 2505.00000
5 2288.16667 2339.66667
6 2525.00000 2288.16667
7 3063.16667 2525.00000
8 2957.83333 3063.16667
9 3222.16667 2957.83333
10 2801.50000 3222.16667
11 2617.00000 2801.50000
12 1347.83333 2617.00000
13 1713.16667 1347.83333
14 1397.16667 1713.16667
15 1924.00000 1397.16667
16 2104.66667 1924.00000
17 1913.16667 2104.66667
18 1661.00000 1913.16667
19 1540.16667 1661.00000
20 2079.83333 1540.16667
21 1206.16667 2079.83333
22 875.50000 1206.16667
23 219.00000 875.50000
24 -152.16667 219.00000
25 -615.83333 -152.16667
26 -991.83333 -615.83333
27 -1167.00000 -991.83333
28 -1048.33333 -1167.00000
29 -1438.83333 -1048.33333
30 -1261.00000 -1438.83333
31 -1387.83333 -1261.00000
32 -1812.16667 -1387.83333
33 -1791.83333 -1812.16667
34 -1819.50000 -1791.83333
35 -1806.00000 -1819.50000
36 -2420.16667 -1806.00000
37 -2625.83333 -2420.16667
38 -2165.83333 -2625.83333
39 -2403.00000 -2165.83333
40 -2556.33333 -2403.00000
41 -2420.83333 -2556.33333
42 -2364.00000 -2420.83333
43 -2494.83333 -2364.00000
44 -2494.16667 -2494.83333
45 -2236.83333 -2494.16667
46 -2026.50000 -2236.83333
47 -1594.00000 -2026.50000
48 -1732.16667 -1594.00000
49 -1312.83333 -1732.16667
50 -920.83333 -1312.83333
51 -888.00000 -920.83333
52 -770.33333 -888.00000
53 -364.83333 -770.33333
54 -474.00000 -364.83333
55 -386.83333 -474.00000
56 -349.16667 -386.83333
57 -73.83333 -349.16667
58 391.50000 -73.83333
59 1002.00000 391.50000
60 93.83333 1002.00000
61 183.16667 93.83333
62 -12.83333 183.16667
63 29.00000 -12.83333
64 -69.33333 29.00000
65 23.16667 -69.33333
66 -87.00000 23.16667
67 -333.83333 -87.00000
68 -382.16667 -333.83333
69 -325.83333 -382.16667
70 -222.50000 -325.83333
71 -438.00000 -222.50000
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/7u98z1322553197.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/8ys2d1322553197.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/9v5qp1322553197.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
hat values (leverages) are all = 0.1666667
and there are no factor predictors; no plot no. 5
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/wessaorg/rcomp/tmp/10pn751322553197.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/wessaorg/rcomp/tmp/11no8c1322553197.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/wessaorg/rcomp/tmp/12tgah1322553197.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/wessaorg/rcomp/tmp/133l0n1322553197.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/wessaorg/rcomp/tmp/14kzee1322553197.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/wessaorg/rcomp/tmp/159me81322553197.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/wessaorg/rcomp/tmp/16fglb1322553197.tab")
+ }
>
> try(system("convert tmp/1e32d1322553197.ps tmp/1e32d1322553197.png",intern=TRUE))
character(0)
> try(system("convert tmp/2tkwi1322553197.ps tmp/2tkwi1322553197.png",intern=TRUE))
character(0)
> try(system("convert tmp/3f0h21322553197.ps tmp/3f0h21322553197.png",intern=TRUE))
character(0)
> try(system("convert tmp/4o48y1322553197.ps tmp/4o48y1322553197.png",intern=TRUE))
character(0)
> try(system("convert tmp/5nnai1322553197.ps tmp/5nnai1322553197.png",intern=TRUE))
character(0)
> try(system("convert tmp/6jzzq1322553197.ps tmp/6jzzq1322553197.png",intern=TRUE))
character(0)
> try(system("convert tmp/7u98z1322553197.ps tmp/7u98z1322553197.png",intern=TRUE))
character(0)
> try(system("convert tmp/8ys2d1322553197.ps tmp/8ys2d1322553197.png",intern=TRUE))
character(0)
> try(system("convert tmp/9v5qp1322553197.ps tmp/9v5qp1322553197.png",intern=TRUE))
character(0)
> try(system("convert tmp/10pn751322553197.ps tmp/10pn751322553197.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.206 0.681 3.951