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 = '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 t
1 16111 1 0 0 0 0 0 0 0 0 0 0 1
2 15554 0 1 0 0 0 0 0 0 0 0 0 2
3 15220 0 0 1 0 0 0 0 0 0 0 0 3
4 14807 0 0 0 1 0 0 0 0 0 0 0 4
5 14291 0 0 0 0 1 0 0 0 0 0 0 5
6 14653 0 0 0 0 0 1 0 0 0 0 0 6
7 17006 0 0 0 0 0 0 1 0 0 0 0 7
8 18032 0 0 0 0 0 0 0 1 0 0 0 8
9 16558 0 0 0 0 0 0 0 0 1 0 0 9
10 16102 0 0 0 0 0 0 0 0 0 1 0 10
11 15055 0 0 0 0 0 0 0 0 0 0 1 11
12 15484 0 0 0 0 0 0 0 0 0 0 0 12
13 14596 1 0 0 0 0 0 0 0 0 0 0 13
14 14609 0 1 0 0 0 0 0 0 0 0 0 14
15 13923 0 0 1 0 0 0 0 0 0 0 0 15
16 14226 0 0 0 1 0 0 0 0 0 0 0 16
17 14056 0 0 0 0 1 0 0 0 0 0 0 17
18 14278 0 0 0 0 0 1 0 0 0 0 0 18
19 16142 0 0 0 0 0 0 1 0 0 0 0 19
20 16509 0 0 0 0 0 0 0 1 0 0 0 20
21 15680 0 0 0 0 0 0 0 0 1 0 0 21
22 14086 0 0 0 0 0 0 0 0 0 1 0 22
23 13129 0 0 0 0 0 0 0 0 0 0 1 23
24 13086 0 0 0 0 0 0 0 0 0 0 0 24
25 13096 1 0 0 0 0 0 0 0 0 0 0 25
26 12280 0 1 0 0 0 0 0 0 0 0 0 26
27 11534 0 0 1 0 0 0 0 0 0 0 0 27
28 11135 0 0 0 1 0 0 0 0 0 0 0 28
29 10903 0 0 0 0 1 0 0 0 0 0 0 29
30 10926 0 0 0 0 0 1 0 0 0 0 0 30
31 13220 0 0 0 0 0 0 1 0 0 0 0 31
32 13581 0 0 0 0 0 0 0 1 0 0 0 32
33 11788 0 0 0 0 0 0 0 0 1 0 0 33
34 11088 0 0 0 0 0 0 0 0 0 1 0 34
35 10434 0 0 0 0 0 0 0 0 0 0 1 35
36 11061 0 0 0 0 0 0 0 0 0 0 0 36
37 10828 1 0 0 0 0 0 0 0 0 0 0 37
38 10270 0 1 0 0 0 0 0 0 0 0 0 38
39 10360 0 0 1 0 0 0 0 0 0 0 0 39
40 9899 0 0 0 1 0 0 0 0 0 0 0 40
41 9395 0 0 0 0 1 0 0 0 0 0 0 41
42 9944 0 0 0 0 0 1 0 0 0 0 0 42
43 12117 0 0 0 0 0 0 1 0 0 0 0 43
44 12474 0 0 0 0 0 0 0 1 0 0 0 44
45 11106 0 0 0 0 0 0 0 0 1 0 0 45
46 10643 0 0 0 0 0 0 0 0 0 1 0 46
47 10227 0 0 0 0 0 0 0 0 0 0 1 47
48 11273 0 0 0 0 0 0 0 0 0 0 0 48
49 11516 1 0 0 0 0 0 0 0 0 0 0 49
50 11583 0 1 0 0 0 0 0 0 0 0 0 50
51 11605 0 0 1 0 0 0 0 0 0 0 0 51
52 11414 0 0 0 1 0 0 0 0 0 0 0 52
53 11181 0 0 0 0 1 0 0 0 0 0 0 53
54 12000 0 0 0 0 0 1 0 0 0 0 0 54
55 14007 0 0 0 0 0 0 1 0 0 0 0 55
56 14582 0 0 0 0 0 0 0 1 0 0 0 56
57 13251 0 0 0 0 0 0 0 0 1 0 0 57
58 12806 0 0 0 0 0 0 0 0 0 1 0 58
59 12645 0 0 0 0 0 0 0 0 0 0 1 59
60 13869 0 0 0 0 0 0 0 0 0 0 0 60
61 13342 1 0 0 0 0 0 0 0 0 0 0 61
62 13079 0 1 0 0 0 0 0 0 0 0 0 62
63 12513 0 0 1 0 0 0 0 0 0 0 0 63
64 12331 0 0 0 1 0 0 0 0 0 0 0 64
65 11882 0 0 0 0 1 0 0 0 0 0 0 65
66 12388 0 0 0 0 0 1 0 0 0 0 0 66
67 14394 0 0 0 0 0 0 1 0 0 0 0 67
68 14635 0 0 0 0 0 0 0 1 0 0 0 68
69 13218 0 0 0 0 0 0 0 0 1 0 0 69
70 12554 0 0 0 0 0 0 0 0 0 1 0 70
71 12031 0 0 0 0 0 0 0 0 0 0 1 71
72 12429 0 0 0 0 0 0 0 0 0 0 0 72
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
14987.65 -174.24 -476.08 -795.59 -968.93 -1269.11
M6 M7 M8 M9 M10 M11
-805.12 1361.54 1899.87 581.69 -88.15 -663.99
t
-50.49
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2323 -1493 533 1180 1911
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 14987.65 729.18 20.554 < 2e-16 ***
M1 -174.24 892.84 -0.195 0.8459
M2 -476.08 891.92 -0.534 0.5955
M3 -795.59 891.09 -0.893 0.3756
M4 -968.93 890.34 -1.088 0.2809
M5 -1269.11 889.68 -1.426 0.1590
M6 -805.12 889.11 -0.906 0.3689
M7 1361.54 888.63 1.532 0.1308
M8 1899.87 888.23 2.139 0.0366 *
M9 581.69 887.92 0.655 0.5149
M10 -88.15 887.70 -0.099 0.9212
M11 -663.99 887.57 -0.748 0.4574
t -50.49 8.84 -5.712 3.87e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1537 on 59 degrees of freedom
Multiple R-squared: 0.488, Adjusted R-squared: 0.3838
F-statistic: 4.686 on 12 and 59 DF, p-value: 2.644e-05
> 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.0116116160 0.023223232 0.98838838
[2,] 0.0098949389 0.019789878 0.99010506
[3,] 0.0044811221 0.008962244 0.99551888
[4,] 0.0014027103 0.002805421 0.99859729
[5,] 0.0011825910 0.002365182 0.99881741
[6,] 0.0007523107 0.001504621 0.99924769
[7,] 0.0031072877 0.006214575 0.99689271
[8,] 0.0063911434 0.012782287 0.99360886
[9,] 0.0200451126 0.040090225 0.97995489
[10,] 0.0230096574 0.046019315 0.97699034
[11,] 0.0380167969 0.076033594 0.96198320
[12,] 0.0560538112 0.112107622 0.94394619
[13,] 0.0985193117 0.197038623 0.90148069
[14,] 0.1424822238 0.284964448 0.85751778
[15,] 0.1865959402 0.373191880 0.81340406
[16,] 0.2275420911 0.455084182 0.77245791
[17,] 0.3305326052 0.661065210 0.66946739
[18,] 0.5154629065 0.969074187 0.48453709
[19,] 0.5948428435 0.810314313 0.40515716
[20,] 0.5979292318 0.804141536 0.40207077
[21,] 0.5539104201 0.892179160 0.44608958
[22,] 0.4724434225 0.944886845 0.52755658
[23,] 0.4021037438 0.804207488 0.59789626
[24,] 0.3570944378 0.714188876 0.64290556
[25,] 0.2982168475 0.596433695 0.70178315
[26,] 0.2485716418 0.497143284 0.75142836
[27,] 0.2268939031 0.453787806 0.77310610
[28,] 0.1978899015 0.395779803 0.80211010
[29,] 0.1743737095 0.348747419 0.82562629
[30,] 0.1620933268 0.324186654 0.83790667
[31,] 0.1695971453 0.339194291 0.83040285
[32,] 0.2489200076 0.497840015 0.75107999
[33,] 0.4333642084 0.866728417 0.56663579
[34,] 0.7396678313 0.520664337 0.26033217
[35,] 0.9078855599 0.184228880 0.09211444
[36,] 0.9474746199 0.105050760 0.05252538
[37,] 0.9677570360 0.064485928 0.03224296
[38,] 0.9740900283 0.051819943 0.02590997
[39,] 0.9703653969 0.059269206 0.02963460
[40,] 0.9652427168 0.069514566 0.03475728
[41,] 0.9368294534 0.126341093 0.06317055
> postscript(file="/var/wessaorg/rcomp/tmp/1ze051322508281.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/2e4n61322508281.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/3wf701322508281.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/4m48l1322508281.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/519fg1322508281.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
1348.08333 1143.41667 1179.41667 990.25000 824.91667 773.41667
7 8 9 10 11 12
1010.25000 1548.41667 1443.08333 1707.41667 1286.75000 1102.25000
13 14 15 16 17 18
438.98333 804.31667 488.31667 1015.15000 1195.81667 1004.31667
19 20 21 22 23 24
752.15000 631.31667 1170.98333 297.31667 -33.35000 -689.85000
25 26 27 28 29 30
-455.11667 -918.78333 -1294.78333 -1469.95000 -1351.28333 -1741.78333
31 32 33 34 35 36
-1563.95000 -1690.78333 -2115.11667 -2094.78333 -2122.45000 -2108.95000
37 38 39 40 41 42
-2117.21667 -2322.88333 -1862.88333 -2100.05000 -2253.38333 -2117.88333
43 44 45 46 47 48
-2061.05000 -2191.88333 -2191.21667 -1933.88333 -1723.55000 -1291.05000
49 50 51 52 53 54
-823.31667 -403.98333 -11.98333 20.85000 138.51667 544.01667
55 56 57 58 59 60
434.85000 522.01667 559.68333 835.01667 1300.35000 1910.85000
61 62 63 64 65 66
1608.58333 1697.91667 1501.91667 1543.75000 1445.41667 1537.91667
67 68 69 70 71 72
1427.75000 1180.91667 1132.58333 1188.91667 1292.25000 1076.75000
> postscript(file="/var/wessaorg/rcomp/tmp/66c4e1322508281.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 1348.08333 NA
1 1143.41667 1348.08333
2 1179.41667 1143.41667
3 990.25000 1179.41667
4 824.91667 990.25000
5 773.41667 824.91667
6 1010.25000 773.41667
7 1548.41667 1010.25000
8 1443.08333 1548.41667
9 1707.41667 1443.08333
10 1286.75000 1707.41667
11 1102.25000 1286.75000
12 438.98333 1102.25000
13 804.31667 438.98333
14 488.31667 804.31667
15 1015.15000 488.31667
16 1195.81667 1015.15000
17 1004.31667 1195.81667
18 752.15000 1004.31667
19 631.31667 752.15000
20 1170.98333 631.31667
21 297.31667 1170.98333
22 -33.35000 297.31667
23 -689.85000 -33.35000
24 -455.11667 -689.85000
25 -918.78333 -455.11667
26 -1294.78333 -918.78333
27 -1469.95000 -1294.78333
28 -1351.28333 -1469.95000
29 -1741.78333 -1351.28333
30 -1563.95000 -1741.78333
31 -1690.78333 -1563.95000
32 -2115.11667 -1690.78333
33 -2094.78333 -2115.11667
34 -2122.45000 -2094.78333
35 -2108.95000 -2122.45000
36 -2117.21667 -2108.95000
37 -2322.88333 -2117.21667
38 -1862.88333 -2322.88333
39 -2100.05000 -1862.88333
40 -2253.38333 -2100.05000
41 -2117.88333 -2253.38333
42 -2061.05000 -2117.88333
43 -2191.88333 -2061.05000
44 -2191.21667 -2191.88333
45 -1933.88333 -2191.21667
46 -1723.55000 -1933.88333
47 -1291.05000 -1723.55000
48 -823.31667 -1291.05000
49 -403.98333 -823.31667
50 -11.98333 -403.98333
51 20.85000 -11.98333
52 138.51667 20.85000
53 544.01667 138.51667
54 434.85000 544.01667
55 522.01667 434.85000
56 559.68333 522.01667
57 835.01667 559.68333
58 1300.35000 835.01667
59 1910.85000 1300.35000
60 1608.58333 1910.85000
61 1697.91667 1608.58333
62 1501.91667 1697.91667
63 1543.75000 1501.91667
64 1445.41667 1543.75000
65 1537.91667 1445.41667
66 1427.75000 1537.91667
67 1180.91667 1427.75000
68 1132.58333 1180.91667
69 1188.91667 1132.58333
70 1292.25000 1188.91667
71 1076.75000 1292.25000
72 NA 1076.75000
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 1143.41667 1348.08333
[2,] 1179.41667 1143.41667
[3,] 990.25000 1179.41667
[4,] 824.91667 990.25000
[5,] 773.41667 824.91667
[6,] 1010.25000 773.41667
[7,] 1548.41667 1010.25000
[8,] 1443.08333 1548.41667
[9,] 1707.41667 1443.08333
[10,] 1286.75000 1707.41667
[11,] 1102.25000 1286.75000
[12,] 438.98333 1102.25000
[13,] 804.31667 438.98333
[14,] 488.31667 804.31667
[15,] 1015.15000 488.31667
[16,] 1195.81667 1015.15000
[17,] 1004.31667 1195.81667
[18,] 752.15000 1004.31667
[19,] 631.31667 752.15000
[20,] 1170.98333 631.31667
[21,] 297.31667 1170.98333
[22,] -33.35000 297.31667
[23,] -689.85000 -33.35000
[24,] -455.11667 -689.85000
[25,] -918.78333 -455.11667
[26,] -1294.78333 -918.78333
[27,] -1469.95000 -1294.78333
[28,] -1351.28333 -1469.95000
[29,] -1741.78333 -1351.28333
[30,] -1563.95000 -1741.78333
[31,] -1690.78333 -1563.95000
[32,] -2115.11667 -1690.78333
[33,] -2094.78333 -2115.11667
[34,] -2122.45000 -2094.78333
[35,] -2108.95000 -2122.45000
[36,] -2117.21667 -2108.95000
[37,] -2322.88333 -2117.21667
[38,] -1862.88333 -2322.88333
[39,] -2100.05000 -1862.88333
[40,] -2253.38333 -2100.05000
[41,] -2117.88333 -2253.38333
[42,] -2061.05000 -2117.88333
[43,] -2191.88333 -2061.05000
[44,] -2191.21667 -2191.88333
[45,] -1933.88333 -2191.21667
[46,] -1723.55000 -1933.88333
[47,] -1291.05000 -1723.55000
[48,] -823.31667 -1291.05000
[49,] -403.98333 -823.31667
[50,] -11.98333 -403.98333
[51,] 20.85000 -11.98333
[52,] 138.51667 20.85000
[53,] 544.01667 138.51667
[54,] 434.85000 544.01667
[55,] 522.01667 434.85000
[56,] 559.68333 522.01667
[57,] 835.01667 559.68333
[58,] 1300.35000 835.01667
[59,] 1910.85000 1300.35000
[60,] 1608.58333 1910.85000
[61,] 1697.91667 1608.58333
[62,] 1501.91667 1697.91667
[63,] 1543.75000 1501.91667
[64,] 1445.41667 1543.75000
[65,] 1537.91667 1445.41667
[66,] 1427.75000 1537.91667
[67,] 1180.91667 1427.75000
[68,] 1132.58333 1180.91667
[69,] 1188.91667 1132.58333
[70,] 1292.25000 1188.91667
[71,] 1076.75000 1292.25000
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 1143.41667 1348.08333
2 1179.41667 1143.41667
3 990.25000 1179.41667
4 824.91667 990.25000
5 773.41667 824.91667
6 1010.25000 773.41667
7 1548.41667 1010.25000
8 1443.08333 1548.41667
9 1707.41667 1443.08333
10 1286.75000 1707.41667
11 1102.25000 1286.75000
12 438.98333 1102.25000
13 804.31667 438.98333
14 488.31667 804.31667
15 1015.15000 488.31667
16 1195.81667 1015.15000
17 1004.31667 1195.81667
18 752.15000 1004.31667
19 631.31667 752.15000
20 1170.98333 631.31667
21 297.31667 1170.98333
22 -33.35000 297.31667
23 -689.85000 -33.35000
24 -455.11667 -689.85000
25 -918.78333 -455.11667
26 -1294.78333 -918.78333
27 -1469.95000 -1294.78333
28 -1351.28333 -1469.95000
29 -1741.78333 -1351.28333
30 -1563.95000 -1741.78333
31 -1690.78333 -1563.95000
32 -2115.11667 -1690.78333
33 -2094.78333 -2115.11667
34 -2122.45000 -2094.78333
35 -2108.95000 -2122.45000
36 -2117.21667 -2108.95000
37 -2322.88333 -2117.21667
38 -1862.88333 -2322.88333
39 -2100.05000 -1862.88333
40 -2253.38333 -2100.05000
41 -2117.88333 -2253.38333
42 -2061.05000 -2117.88333
43 -2191.88333 -2061.05000
44 -2191.21667 -2191.88333
45 -1933.88333 -2191.21667
46 -1723.55000 -1933.88333
47 -1291.05000 -1723.55000
48 -823.31667 -1291.05000
49 -403.98333 -823.31667
50 -11.98333 -403.98333
51 20.85000 -11.98333
52 138.51667 20.85000
53 544.01667 138.51667
54 434.85000 544.01667
55 522.01667 434.85000
56 559.68333 522.01667
57 835.01667 559.68333
58 1300.35000 835.01667
59 1910.85000 1300.35000
60 1608.58333 1910.85000
61 1697.91667 1608.58333
62 1501.91667 1697.91667
63 1543.75000 1501.91667
64 1445.41667 1543.75000
65 1537.91667 1445.41667
66 1427.75000 1537.91667
67 1180.91667 1427.75000
68 1132.58333 1180.91667
69 1188.91667 1132.58333
70 1292.25000 1188.91667
71 1076.75000 1292.25000
> 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/7dtbu1322508281.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/8s6ir1322508281.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/9kixy1322508281.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')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/wessaorg/rcomp/tmp/10nk4j1322508281.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/110cjm1322508281.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/12gs7i1322508282.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/131nz81322508282.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/14z3md1322508282.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/15xy7q1322508282.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/164qhp1322508282.tab")
+ }
>
> try(system("convert tmp/1ze051322508281.ps tmp/1ze051322508281.png",intern=TRUE))
character(0)
> try(system("convert tmp/2e4n61322508281.ps tmp/2e4n61322508281.png",intern=TRUE))
character(0)
> try(system("convert tmp/3wf701322508281.ps tmp/3wf701322508281.png",intern=TRUE))
character(0)
> try(system("convert tmp/4m48l1322508281.ps tmp/4m48l1322508281.png",intern=TRUE))
character(0)
> try(system("convert tmp/519fg1322508281.ps tmp/519fg1322508281.png",intern=TRUE))
character(0)
> try(system("convert tmp/66c4e1322508281.ps tmp/66c4e1322508281.png",intern=TRUE))
character(0)
> try(system("convert tmp/7dtbu1322508281.ps tmp/7dtbu1322508281.png",intern=TRUE))
character(0)
> try(system("convert tmp/8s6ir1322508281.ps tmp/8s6ir1322508281.png",intern=TRUE))
character(0)
> try(system("convert tmp/9kixy1322508281.ps tmp/9kixy1322508281.png",intern=TRUE))
character(0)
> try(system("convert tmp/10nk4j1322508281.ps tmp/10nk4j1322508281.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.419 0.511 3.941