R version 2.12.0 (2010-10-15)
Copyright (C) 2010 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(1
+ ,216.234
+ ,627
+ ,1.59
+ ,2
+ ,213.586
+ ,696
+ ,1.26
+ ,3
+ ,209.465
+ ,825
+ ,1.13
+ ,4
+ ,204.045
+ ,677
+ ,1.92
+ ,5
+ ,200.237
+ ,656
+ ,2.61
+ ,6
+ ,203.666
+ ,785
+ ,2.26
+ ,7
+ ,241.476
+ ,412
+ ,2.41
+ ,8
+ ,260.307
+ ,352
+ ,2.26
+ ,9
+ ,243.324
+ ,839
+ ,2.03
+ ,10
+ ,244.460
+ ,729
+ ,2.86
+ ,11
+ ,233.575
+ ,696
+ ,2.55
+ ,12
+ ,237.217
+ ,641
+ ,2.27
+ ,1
+ ,235.243
+ ,695
+ ,2.26
+ ,2
+ ,230.354
+ ,638
+ ,2.57
+ ,3
+ ,227.184
+ ,762
+ ,3.07
+ ,4
+ ,221.678
+ ,635
+ ,2.76
+ ,5
+ ,217.142
+ ,721
+ ,2.51
+ ,6
+ ,219.452
+ ,854
+ ,2.87
+ ,7
+ ,256.446
+ ,418
+ ,3.14
+ ,8
+ ,265.845
+ ,367
+ ,3.11
+ ,9
+ ,248.624
+ ,824
+ ,3.16
+ ,10
+ ,241.114
+ ,687
+ ,2.47
+ ,11
+ ,229.245
+ ,601
+ ,2.57
+ ,12
+ ,231.805
+ ,676
+ ,2.89
+ ,1
+ ,219.277
+ ,740
+ ,2.63
+ ,2
+ ,219.313
+ ,691
+ ,2.38
+ ,3
+ ,212.610
+ ,683
+ ,1.69
+ ,4
+ ,214.771
+ ,594
+ ,1.96
+ ,5
+ ,211.142
+ ,729
+ ,2.19
+ ,6
+ ,211.457
+ ,731
+ ,1.87
+ ,7
+ ,240.048
+ ,386
+ ,1.6
+ ,8
+ ,240.636
+ ,331
+ ,1.63
+ ,9
+ ,230.580
+ ,707
+ ,1.22
+ ,10
+ ,208.795
+ ,715
+ ,1.21
+ ,11
+ ,197.922
+ ,657
+ ,1.49
+ ,12
+ ,194.596
+ ,653
+ ,1.64
+ ,1
+ ,194.581
+ ,642
+ ,1.66
+ ,2
+ ,185.686
+ ,643
+ ,1.77
+ ,3
+ ,178.106
+ ,718
+ ,1.82
+ ,4
+ ,172.608
+ ,654
+ ,1.78
+ ,5
+ ,167.302
+ ,632
+ ,1.28
+ ,6
+ ,168.053
+ ,731
+ ,1.29
+ ,7
+ ,202.300
+ ,392
+ ,1.37
+ ,8
+ ,202.388
+ ,344
+ ,1.12
+ ,9
+ ,182.516
+ ,792
+ ,1.51
+ ,10
+ ,173.476
+ ,852
+ ,2.24
+ ,11
+ ,166.444
+ ,649
+ ,2.94
+ ,12
+ ,171.297
+ ,629
+ ,3.09
+ ,1
+ ,169.701
+ ,685
+ ,3.46
+ ,2
+ ,164.182
+ ,617
+ ,3.64
+ ,3
+ ,161.914
+ ,715
+ ,4.39
+ ,4
+ ,159.612
+ ,715
+ ,4.15
+ ,5
+ ,151.001
+ ,629
+ ,5.21
+ ,6
+ ,158.114
+ ,916
+ ,5.8
+ ,7
+ ,186.530
+ ,531
+ ,5.91
+ ,8
+ ,187.069
+ ,357
+ ,5.39
+ ,9
+ ,174.330
+ ,917
+ ,5.46
+ ,10
+ ,169.362
+ ,828
+ ,4.72
+ ,11
+ ,166.827
+ ,708
+ ,3.14
+ ,12
+ ,178.037
+ ,858
+ ,2.63
+ ,1
+ ,186.413
+ ,775
+ ,2.32
+ ,2
+ ,189.226
+ ,785
+ ,1.93
+ ,3
+ ,191.563
+ ,1006
+ ,0.62
+ ,4
+ ,188.906
+ ,789
+ ,0.6
+ ,5
+ ,186.005
+ ,734
+ ,-0.37
+ ,6
+ ,195.309
+ ,906
+ ,-1.1
+ ,7
+ ,223.532
+ ,532
+ ,-1.68
+ ,8
+ ,226.899
+ ,387
+ ,-0.78
+ ,9
+ ,214.126
+ ,991
+ ,-1.19
+ ,10
+ ,206.903
+ ,841
+ ,-0.97
+ ,11
+ ,204.442
+ ,892
+ ,-0.12
+ ,12
+ ,220.375
+ ,782
+ ,0.26)
+ ,dim=c(4
+ ,72)
+ ,dimnames=list(c('month'
+ ,'werklozen'
+ ,'faillissementen'
+ ,'inflatie')
+ ,1:72))
> y <- array(NA,dim=c(4,72),dimnames=list(c('month','werklozen','faillissementen','inflatie'),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 = 'Do not include Seasonal Dummies'
> par1 = '2'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> 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
werklozen month faillissementen inflatie
1 216.234 1 627 1.59
2 213.586 2 696 1.26
3 209.465 3 825 1.13
4 204.045 4 677 1.92
5 200.237 5 656 2.61
6 203.666 6 785 2.26
7 241.476 7 412 2.41
8 260.307 8 352 2.26
9 243.324 9 839 2.03
10 244.460 10 729 2.86
11 233.575 11 696 2.55
12 237.217 12 641 2.27
13 235.243 1 695 2.26
14 230.354 2 638 2.57
15 227.184 3 762 3.07
16 221.678 4 635 2.76
17 217.142 5 721 2.51
18 219.452 6 854 2.87
19 256.446 7 418 3.14
20 265.845 8 367 3.11
21 248.624 9 824 3.16
22 241.114 10 687 2.47
23 229.245 11 601 2.57
24 231.805 12 676 2.89
25 219.277 1 740 2.63
26 219.313 2 691 2.38
27 212.610 3 683 1.69
28 214.771 4 594 1.96
29 211.142 5 729 2.19
30 211.457 6 731 1.87
31 240.048 7 386 1.60
32 240.636 8 331 1.63
33 230.580 9 707 1.22
34 208.795 10 715 1.21
35 197.922 11 657 1.49
36 194.596 12 653 1.64
37 194.581 1 642 1.66
38 185.686 2 643 1.77
39 178.106 3 718 1.82
40 172.608 4 654 1.78
41 167.302 5 632 1.28
42 168.053 6 731 1.29
43 202.300 7 392 1.37
44 202.388 8 344 1.12
45 182.516 9 792 1.51
46 173.476 10 852 2.24
47 166.444 11 649 2.94
48 171.297 12 629 3.09
49 169.701 1 685 3.46
50 164.182 2 617 3.64
51 161.914 3 715 4.39
52 159.612 4 715 4.15
53 151.001 5 629 5.21
54 158.114 6 916 5.80
55 186.530 7 531 5.91
56 187.069 8 357 5.39
57 174.330 9 917 5.46
58 169.362 10 828 4.72
59 166.827 11 708 3.14
60 178.037 12 858 2.63
61 186.413 1 775 2.32
62 189.226 2 785 1.93
63 191.563 3 1006 0.62
64 188.906 4 789 0.60
65 186.005 5 734 -0.37
66 195.309 6 906 -1.10
67 223.532 7 532 -1.68
68 226.899 8 387 -0.78
69 214.126 9 991 -1.19
70 206.903 10 841 -0.97
71 204.442 11 892 -0.12
72 220.375 12 782 0.26
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) month faillissementen inflatie
252.40736 1.10477 -0.06416 -5.07414
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-43.58 -20.98 -1.92 19.64 55.18
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 252.40736 15.73979 16.036 < 2e-16 ***
month 1.10477 0.86887 1.272 0.20788
faillissementen -0.06416 0.01936 -3.314 0.00148 **
inflatie -5.07414 1.94836 -2.604 0.01130 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 25.42 on 68 degrees of freedom
Multiple R-squared: 0.2128, Adjusted R-squared: 0.1781
F-statistic: 6.128 on 3 and 68 DF, p-value: 0.0009435
> 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.012126796 0.0242535927 9.878732e-01
[2,] 0.002830696 0.0056613928 9.971693e-01
[3,] 0.011776950 0.0235538994 9.882231e-01
[4,] 0.011729080 0.0234581599 9.882709e-01
[5,] 0.012524844 0.0250496878 9.874752e-01
[6,] 0.019971801 0.0399436030 9.800282e-01
[7,] 0.078093013 0.1561860256 9.219070e-01
[8,] 0.059756990 0.1195139791 9.402430e-01
[9,] 0.044306110 0.0886122204 9.556939e-01
[10,] 0.028694257 0.0573885149 9.713057e-01
[11,] 0.018268136 0.0365362720 9.817319e-01
[12,] 0.012162185 0.0243243694 9.878378e-01
[13,] 0.011335403 0.0226708058 9.886646e-01
[14,] 0.014918247 0.0298364942 9.850818e-01
[15,] 0.037742086 0.0754841729 9.622579e-01
[16,] 0.042621977 0.0852439539 9.573780e-01
[17,] 0.053955565 0.1079111297 9.460444e-01
[18,] 0.076216263 0.1524325267 9.237837e-01
[19,] 0.083101811 0.1662036223 9.168982e-01
[20,] 0.090822540 0.1816450790 9.091775e-01
[21,] 0.080724795 0.1614495902 9.192752e-01
[22,] 0.081974826 0.1639496519 9.180252e-01
[23,] 0.092735602 0.1854712046 9.072644e-01
[24,] 0.096341297 0.1926825932 9.036587e-01
[25,] 0.136127392 0.2722547830 8.638726e-01
[26,] 0.224360125 0.4487202496 7.756399e-01
[27,] 0.418430839 0.8368616781 5.815692e-01
[28,] 0.438993828 0.8779876556 5.610062e-01
[29,] 0.562537889 0.8749242225 4.374621e-01
[30,] 0.676805986 0.6463880285 3.231940e-01
[31,] 0.698537024 0.6029259520 3.014630e-01
[32,] 0.752457446 0.4950851078 2.475426e-01
[33,] 0.809590208 0.3808195835 1.904098e-01
[34,] 0.900480093 0.1990398139 9.951991e-02
[35,] 0.964906551 0.0701868973 3.509345e-02
[36,] 0.984539857 0.0309202858 1.546014e-02
[37,] 0.978891336 0.0422173279 2.110866e-02
[38,] 0.969688016 0.0606239677 3.031198e-02
[39,] 0.963031492 0.0739370163 3.696851e-02
[40,] 0.975288474 0.0494230527 2.471153e-02
[41,] 0.997503980 0.0049920397 2.496020e-03
[42,] 0.999535695 0.0009286090 4.643045e-04
[43,] 0.999638785 0.0007224306 3.612153e-04
[44,] 0.999758032 0.0004839362 2.419681e-04
[45,] 0.999726014 0.0005479729 2.739864e-04
[46,] 0.999731849 0.0005363026 2.681513e-04
[47,] 0.999904062 0.0001918750 9.593752e-05
[48,] 0.999759622 0.0004807551 2.403776e-04
[49,] 0.999621269 0.0007574613 3.787306e-04
[50,] 0.999198895 0.0016022102 8.011051e-04
[51,] 0.998998642 0.0020027160 1.001358e-03
[52,] 0.997686175 0.0046276503 2.313825e-03
[53,] 0.998230247 0.0035395052 1.769753e-03
[54,] 0.999094306 0.0018113883 9.056942e-04
[55,] 0.997153256 0.0056934889 2.846744e-03
[56,] 0.992317737 0.0153645256 7.682263e-03
[57,] 0.993373494 0.0132530125 6.626506e-03
[58,] 0.990038086 0.0199238282 9.961914e-03
[59,] 0.966412977 0.0671740468 3.358702e-02
> postscript(file="/var/www/rcomp/tmp/1eeqf1292944729.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/www/rcomp/tmp/2eeqf1292944729.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/www/rcomp/tmp/37nqi1292944729.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/www/rcomp/tmp/47nqi1292944729.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/www/rcomp/tmp/57nqi1292944729.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
11.0192811 10.0192147 12.4106925 0.3985264 -2.3604888 6.4646778
7 8 9 10 11 12
19.9986263 33.1140164 45.1060512 42.2910041 26.6109036 24.1984649
13 14 15 16 17 18
37.7909676 29.7129489 35.9313266 19.5990046 18.2076224 29.7730777
19 20 21 22 23 24
39.0577216 43.9274662 55.1774028 34.2712876 16.2870030 24.1781006
25 26 27 28 29 30
26.5896870 21.1084442 9.2862180 6.0020518 11.0971924 8.8120181
31 32 33 34 35 36
12.7923609 8.8989063 19.7826212 -2.6445975 -16.9230024 -20.8493015
37 38 39 40 41 42
-9.3161000 -18.6935552 -22.3124760 -33.2245779 -43.5839842 -37.5349844
43 44 45 46 47 48
-25.7377202 -31.1028011 -21.3561135 -23.9470466 -41.5567917 -38.3306817
49 50 51 52 53 54
-22.3036812 -32.3771196 -25.6564165 -30.2809834 -40.1360915 -12.7196462
55 56 57 58 59 60
-9.5526066 -23.9211091 -1.4790108 -17.0170605 -36.3734094 -19.2317050
61 62 63 64 65 66
-5.6016295 -5.2306983 3.5341895 -14.2522050 -26.7088018 -11.1778462
67 68 69 70 71 72
-10.9991833 -13.4737081 9.3219270 -7.5138240 -3.4943175 6.2042714
> postscript(file="/var/www/rcomp/tmp/60xp31292944729.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 11.0192811 NA
1 10.0192147 11.0192811
2 12.4106925 10.0192147
3 0.3985264 12.4106925
4 -2.3604888 0.3985264
5 6.4646778 -2.3604888
6 19.9986263 6.4646778
7 33.1140164 19.9986263
8 45.1060512 33.1140164
9 42.2910041 45.1060512
10 26.6109036 42.2910041
11 24.1984649 26.6109036
12 37.7909676 24.1984649
13 29.7129489 37.7909676
14 35.9313266 29.7129489
15 19.5990046 35.9313266
16 18.2076224 19.5990046
17 29.7730777 18.2076224
18 39.0577216 29.7730777
19 43.9274662 39.0577216
20 55.1774028 43.9274662
21 34.2712876 55.1774028
22 16.2870030 34.2712876
23 24.1781006 16.2870030
24 26.5896870 24.1781006
25 21.1084442 26.5896870
26 9.2862180 21.1084442
27 6.0020518 9.2862180
28 11.0971924 6.0020518
29 8.8120181 11.0971924
30 12.7923609 8.8120181
31 8.8989063 12.7923609
32 19.7826212 8.8989063
33 -2.6445975 19.7826212
34 -16.9230024 -2.6445975
35 -20.8493015 -16.9230024
36 -9.3161000 -20.8493015
37 -18.6935552 -9.3161000
38 -22.3124760 -18.6935552
39 -33.2245779 -22.3124760
40 -43.5839842 -33.2245779
41 -37.5349844 -43.5839842
42 -25.7377202 -37.5349844
43 -31.1028011 -25.7377202
44 -21.3561135 -31.1028011
45 -23.9470466 -21.3561135
46 -41.5567917 -23.9470466
47 -38.3306817 -41.5567917
48 -22.3036812 -38.3306817
49 -32.3771196 -22.3036812
50 -25.6564165 -32.3771196
51 -30.2809834 -25.6564165
52 -40.1360915 -30.2809834
53 -12.7196462 -40.1360915
54 -9.5526066 -12.7196462
55 -23.9211091 -9.5526066
56 -1.4790108 -23.9211091
57 -17.0170605 -1.4790108
58 -36.3734094 -17.0170605
59 -19.2317050 -36.3734094
60 -5.6016295 -19.2317050
61 -5.2306983 -5.6016295
62 3.5341895 -5.2306983
63 -14.2522050 3.5341895
64 -26.7088018 -14.2522050
65 -11.1778462 -26.7088018
66 -10.9991833 -11.1778462
67 -13.4737081 -10.9991833
68 9.3219270 -13.4737081
69 -7.5138240 9.3219270
70 -3.4943175 -7.5138240
71 6.2042714 -3.4943175
72 NA 6.2042714
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 10.0192147 11.0192811
[2,] 12.4106925 10.0192147
[3,] 0.3985264 12.4106925
[4,] -2.3604888 0.3985264
[5,] 6.4646778 -2.3604888
[6,] 19.9986263 6.4646778
[7,] 33.1140164 19.9986263
[8,] 45.1060512 33.1140164
[9,] 42.2910041 45.1060512
[10,] 26.6109036 42.2910041
[11,] 24.1984649 26.6109036
[12,] 37.7909676 24.1984649
[13,] 29.7129489 37.7909676
[14,] 35.9313266 29.7129489
[15,] 19.5990046 35.9313266
[16,] 18.2076224 19.5990046
[17,] 29.7730777 18.2076224
[18,] 39.0577216 29.7730777
[19,] 43.9274662 39.0577216
[20,] 55.1774028 43.9274662
[21,] 34.2712876 55.1774028
[22,] 16.2870030 34.2712876
[23,] 24.1781006 16.2870030
[24,] 26.5896870 24.1781006
[25,] 21.1084442 26.5896870
[26,] 9.2862180 21.1084442
[27,] 6.0020518 9.2862180
[28,] 11.0971924 6.0020518
[29,] 8.8120181 11.0971924
[30,] 12.7923609 8.8120181
[31,] 8.8989063 12.7923609
[32,] 19.7826212 8.8989063
[33,] -2.6445975 19.7826212
[34,] -16.9230024 -2.6445975
[35,] -20.8493015 -16.9230024
[36,] -9.3161000 -20.8493015
[37,] -18.6935552 -9.3161000
[38,] -22.3124760 -18.6935552
[39,] -33.2245779 -22.3124760
[40,] -43.5839842 -33.2245779
[41,] -37.5349844 -43.5839842
[42,] -25.7377202 -37.5349844
[43,] -31.1028011 -25.7377202
[44,] -21.3561135 -31.1028011
[45,] -23.9470466 -21.3561135
[46,] -41.5567917 -23.9470466
[47,] -38.3306817 -41.5567917
[48,] -22.3036812 -38.3306817
[49,] -32.3771196 -22.3036812
[50,] -25.6564165 -32.3771196
[51,] -30.2809834 -25.6564165
[52,] -40.1360915 -30.2809834
[53,] -12.7196462 -40.1360915
[54,] -9.5526066 -12.7196462
[55,] -23.9211091 -9.5526066
[56,] -1.4790108 -23.9211091
[57,] -17.0170605 -1.4790108
[58,] -36.3734094 -17.0170605
[59,] -19.2317050 -36.3734094
[60,] -5.6016295 -19.2317050
[61,] -5.2306983 -5.6016295
[62,] 3.5341895 -5.2306983
[63,] -14.2522050 3.5341895
[64,] -26.7088018 -14.2522050
[65,] -11.1778462 -26.7088018
[66,] -10.9991833 -11.1778462
[67,] -13.4737081 -10.9991833
[68,] 9.3219270 -13.4737081
[69,] -7.5138240 9.3219270
[70,] -3.4943175 -7.5138240
[71,] 6.2042714 -3.4943175
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 10.0192147 11.0192811
2 12.4106925 10.0192147
3 0.3985264 12.4106925
4 -2.3604888 0.3985264
5 6.4646778 -2.3604888
6 19.9986263 6.4646778
7 33.1140164 19.9986263
8 45.1060512 33.1140164
9 42.2910041 45.1060512
10 26.6109036 42.2910041
11 24.1984649 26.6109036
12 37.7909676 24.1984649
13 29.7129489 37.7909676
14 35.9313266 29.7129489
15 19.5990046 35.9313266
16 18.2076224 19.5990046
17 29.7730777 18.2076224
18 39.0577216 29.7730777
19 43.9274662 39.0577216
20 55.1774028 43.9274662
21 34.2712876 55.1774028
22 16.2870030 34.2712876
23 24.1781006 16.2870030
24 26.5896870 24.1781006
25 21.1084442 26.5896870
26 9.2862180 21.1084442
27 6.0020518 9.2862180
28 11.0971924 6.0020518
29 8.8120181 11.0971924
30 12.7923609 8.8120181
31 8.8989063 12.7923609
32 19.7826212 8.8989063
33 -2.6445975 19.7826212
34 -16.9230024 -2.6445975
35 -20.8493015 -16.9230024
36 -9.3161000 -20.8493015
37 -18.6935552 -9.3161000
38 -22.3124760 -18.6935552
39 -33.2245779 -22.3124760
40 -43.5839842 -33.2245779
41 -37.5349844 -43.5839842
42 -25.7377202 -37.5349844
43 -31.1028011 -25.7377202
44 -21.3561135 -31.1028011
45 -23.9470466 -21.3561135
46 -41.5567917 -23.9470466
47 -38.3306817 -41.5567917
48 -22.3036812 -38.3306817
49 -32.3771196 -22.3036812
50 -25.6564165 -32.3771196
51 -30.2809834 -25.6564165
52 -40.1360915 -30.2809834
53 -12.7196462 -40.1360915
54 -9.5526066 -12.7196462
55 -23.9211091 -9.5526066
56 -1.4790108 -23.9211091
57 -17.0170605 -1.4790108
58 -36.3734094 -17.0170605
59 -19.2317050 -36.3734094
60 -5.6016295 -19.2317050
61 -5.2306983 -5.6016295
62 3.5341895 -5.2306983
63 -14.2522050 3.5341895
64 -26.7088018 -14.2522050
65 -11.1778462 -26.7088018
66 -10.9991833 -11.1778462
67 -13.4737081 -10.9991833
68 9.3219270 -13.4737081
69 -7.5138240 9.3219270
70 -3.4943175 -7.5138240
71 6.2042714 -3.4943175
> 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/www/rcomp/tmp/7s6o61292944729.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/www/rcomp/tmp/8s6o61292944729.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/www/rcomp/tmp/9s6o61292944729.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/www/rcomp/tmp/10lxor1292944729.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/www/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/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/www/rcomp/tmp/11oy4f1292944729.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/www/rcomp/tmp/12syk21292944729.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/www/rcomp/tmp/13oqib1292944729.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/www/rcomp/tmp/149qzh1292944729.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/www/rcomp/tmp/15vrx51292944729.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/www/rcomp/tmp/16g9wb1292944729.tab")
+ }
>
> try(system("convert tmp/1eeqf1292944729.ps tmp/1eeqf1292944729.png",intern=TRUE))
character(0)
> try(system("convert tmp/2eeqf1292944729.ps tmp/2eeqf1292944729.png",intern=TRUE))
character(0)
> try(system("convert tmp/37nqi1292944729.ps tmp/37nqi1292944729.png",intern=TRUE))
character(0)
> try(system("convert tmp/47nqi1292944729.ps tmp/47nqi1292944729.png",intern=TRUE))
character(0)
> try(system("convert tmp/57nqi1292944729.ps tmp/57nqi1292944729.png",intern=TRUE))
character(0)
> try(system("convert tmp/60xp31292944729.ps tmp/60xp31292944729.png",intern=TRUE))
character(0)
> try(system("convert tmp/7s6o61292944729.ps tmp/7s6o61292944729.png",intern=TRUE))
character(0)
> try(system("convert tmp/8s6o61292944729.ps tmp/8s6o61292944729.png",intern=TRUE))
character(0)
> try(system("convert tmp/9s6o61292944729.ps tmp/9s6o61292944729.png",intern=TRUE))
character(0)
> try(system("convert tmp/10lxor1292944729.ps tmp/10lxor1292944729.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.170 1.560 4.834