R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
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(474605,0,470390,0,461251,0,454724,0,455626,0,516847,0,525192,0,522975,0,518585,0,509239,0,512238,0,519164,0,517009,0,509933,0,509127,0,500875,0,506971,0,569323,0,579714,0,577992,0,565644,0,547344,0,554788,0,562325,0,560854,0,555332,0,543599,0,536662,0,542722,0,593530,0,610763,0,612613,0,611324,0,594167,0,595454,0,590865,0,589379,0,584428,0,573100,0,567456,0,569028,0,620735,0,628884,0,628232,0,612117,0,595404,0,597141,0,593408,0,590072,0,579799,0,574205,0,572775,0,572942,0,619567,0,625809,0,619916,0,587625,0,565724,0,557274,0,560576,0,548854,0,531673,0,525919,0,511038,0,498662,0,555362,0,564591,0,541667,0,527070,0,509846,0,514258,0,516922,0,507561,0,492622,0,490243,0,469357,0,477580,0,528379,0,533590,0,517945,1,506174,1,501866,1,516441,1,528222,1,532638,1),dim=c(2,85),dimnames=list(c('Werkzoekend','Crisis'),1:85))
> y <- array(NA,dim=c(2,85),dimnames=list(c('Werkzoekend','Crisis'),1:85))
> 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'
> #'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
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> 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
Werkzoekend Crisis M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 474605 0 1 0 0 0 0 0 0 0 0 0 0
2 470390 0 0 1 0 0 0 0 0 0 0 0 0
3 461251 0 0 0 1 0 0 0 0 0 0 0 0
4 454724 0 0 0 0 1 0 0 0 0 0 0 0
5 455626 0 0 0 0 0 1 0 0 0 0 0 0
6 516847 0 0 0 0 0 0 1 0 0 0 0 0
7 525192 0 0 0 0 0 0 0 1 0 0 0 0
8 522975 0 0 0 0 0 0 0 0 1 0 0 0
9 518585 0 0 0 0 0 0 0 0 0 1 0 0
10 509239 0 0 0 0 0 0 0 0 0 0 1 0
11 512238 0 0 0 0 0 0 0 0 0 0 0 1
12 519164 0 0 0 0 0 0 0 0 0 0 0 0
13 517009 0 1 0 0 0 0 0 0 0 0 0 0
14 509933 0 0 1 0 0 0 0 0 0 0 0 0
15 509127 0 0 0 1 0 0 0 0 0 0 0 0
16 500875 0 0 0 0 1 0 0 0 0 0 0 0
17 506971 0 0 0 0 0 1 0 0 0 0 0 0
18 569323 0 0 0 0 0 0 1 0 0 0 0 0
19 579714 0 0 0 0 0 0 0 1 0 0 0 0
20 577992 0 0 0 0 0 0 0 0 1 0 0 0
21 565644 0 0 0 0 0 0 0 0 0 1 0 0
22 547344 0 0 0 0 0 0 0 0 0 0 1 0
23 554788 0 0 0 0 0 0 0 0 0 0 0 1
24 562325 0 0 0 0 0 0 0 0 0 0 0 0
25 560854 0 1 0 0 0 0 0 0 0 0 0 0
26 555332 0 0 1 0 0 0 0 0 0 0 0 0
27 543599 0 0 0 1 0 0 0 0 0 0 0 0
28 536662 0 0 0 0 1 0 0 0 0 0 0 0
29 542722 0 0 0 0 0 1 0 0 0 0 0 0
30 593530 0 0 0 0 0 0 1 0 0 0 0 0
31 610763 0 0 0 0 0 0 0 1 0 0 0 0
32 612613 0 0 0 0 0 0 0 0 1 0 0 0
33 611324 0 0 0 0 0 0 0 0 0 1 0 0
34 594167 0 0 0 0 0 0 0 0 0 0 1 0
35 595454 0 0 0 0 0 0 0 0 0 0 0 1
36 590865 0 0 0 0 0 0 0 0 0 0 0 0
37 589379 0 1 0 0 0 0 0 0 0 0 0 0
38 584428 0 0 1 0 0 0 0 0 0 0 0 0
39 573100 0 0 0 1 0 0 0 0 0 0 0 0
40 567456 0 0 0 0 1 0 0 0 0 0 0 0
41 569028 0 0 0 0 0 1 0 0 0 0 0 0
42 620735 0 0 0 0 0 0 1 0 0 0 0 0
43 628884 0 0 0 0 0 0 0 1 0 0 0 0
44 628232 0 0 0 0 0 0 0 0 1 0 0 0
45 612117 0 0 0 0 0 0 0 0 0 1 0 0
46 595404 0 0 0 0 0 0 0 0 0 0 1 0
47 597141 0 0 0 0 0 0 0 0 0 0 0 1
48 593408 0 0 0 0 0 0 0 0 0 0 0 0
49 590072 0 1 0 0 0 0 0 0 0 0 0 0
50 579799 0 0 1 0 0 0 0 0 0 0 0 0
51 574205 0 0 0 1 0 0 0 0 0 0 0 0
52 572775 0 0 0 0 1 0 0 0 0 0 0 0
53 572942 0 0 0 0 0 1 0 0 0 0 0 0
54 619567 0 0 0 0 0 0 1 0 0 0 0 0
55 625809 0 0 0 0 0 0 0 1 0 0 0 0
56 619916 0 0 0 0 0 0 0 0 1 0 0 0
57 587625 0 0 0 0 0 0 0 0 0 1 0 0
58 565724 0 0 0 0 0 0 0 0 0 0 1 0
59 557274 0 0 0 0 0 0 0 0 0 0 0 1
60 560576 0 0 0 0 0 0 0 0 0 0 0 0
61 548854 0 1 0 0 0 0 0 0 0 0 0 0
62 531673 0 0 1 0 0 0 0 0 0 0 0 0
63 525919 0 0 0 1 0 0 0 0 0 0 0 0
64 511038 0 0 0 0 1 0 0 0 0 0 0 0
65 498662 0 0 0 0 0 1 0 0 0 0 0 0
66 555362 0 0 0 0 0 0 1 0 0 0 0 0
67 564591 0 0 0 0 0 0 0 1 0 0 0 0
68 541667 0 0 0 0 0 0 0 0 1 0 0 0
69 527070 0 0 0 0 0 0 0 0 0 1 0 0
70 509846 0 0 0 0 0 0 0 0 0 0 1 0
71 514258 0 0 0 0 0 0 0 0 0 0 0 1
72 516922 0 0 0 0 0 0 0 0 0 0 0 0
73 507561 0 1 0 0 0 0 0 0 0 0 0 0
74 492622 0 0 1 0 0 0 0 0 0 0 0 0
75 490243 0 0 0 1 0 0 0 0 0 0 0 0
76 469357 0 0 0 0 1 0 0 0 0 0 0 0
77 477580 0 0 0 0 0 1 0 0 0 0 0 0
78 528379 0 0 0 0 0 0 1 0 0 0 0 0
79 533590 0 0 0 0 0 0 0 1 0 0 0 0
80 517945 1 0 0 0 0 0 0 0 1 0 0 0
81 506174 1 0 0 0 0 0 0 0 0 1 0 0
82 501866 1 0 0 0 0 0 0 0 0 0 1 0
83 516441 1 0 0 0 0 0 0 0 0 0 0 1
84 528222 1 0 0 0 0 0 0 0 0 0 0 0
85 532638 1 1 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) Crisis M1 M2 M3 M4
559200 -42917 -13714 -27175 -33851 -43073
M5 M6 M7 M8 M9 M10
-41553 12763 22021 21408 8151 -6842
M11
-3413
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-70881.2 -38941.2 -352.3 39307.8 56648.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 559200 15588 35.874 <2e-16 ***
Crisis -42918 17908 -2.397 0.0191 *
M1 -13714 21058 -0.651 0.5170
M2 -27175 21896 -1.241 0.2186
M3 -33851 21896 -1.546 0.1265
M4 -43073 21896 -1.967 0.0530 .
M5 -41553 21896 -1.898 0.0617 .
M6 12763 21896 0.583 0.5618
M7 22020 21896 1.006 0.3179
M8 21408 21746 0.984 0.3282
M9 8151 21746 0.375 0.7089
M10 -6842 21746 -0.315 0.7540
M11 -3413 21746 -0.157 0.8757
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 40680 on 72 degrees of freedom
Multiple R-squared: 0.2829, Adjusted R-squared: 0.1633
F-statistic: 2.366 on 12 and 72 DF, p-value: 0.01250
> 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.6174032 0.765193587 0.382596793
[2,] 0.6009651 0.798069734 0.399034867
[3,] 0.5890556 0.821888871 0.410944435
[4,] 0.5848507 0.830298697 0.415149349
[5,] 0.5809523 0.838095450 0.419047725
[6,] 0.5468054 0.906389128 0.453194564
[7,] 0.4896617 0.979323362 0.510338319
[8,] 0.4457627 0.891525425 0.554237288
[9,] 0.4045312 0.809062325 0.595468838
[10,] 0.4671451 0.934290103 0.532854948
[11,] 0.5160494 0.967901185 0.483950593
[12,] 0.5287658 0.942468362 0.471234181
[13,] 0.5377884 0.924423166 0.462211583
[14,] 0.5509860 0.898027956 0.449013978
[15,] 0.5306298 0.938740333 0.469370166
[16,] 0.5319174 0.936165119 0.468082559
[17,] 0.5392084 0.921583248 0.460791624
[18,] 0.5735167 0.852966650 0.426483325
[19,] 0.5929774 0.814045107 0.407022554
[20,] 0.5993530 0.801294009 0.400647005
[21,] 0.5759896 0.848020895 0.424010447
[22,] 0.6049121 0.790175823 0.395087912
[23,] 0.6458326 0.708334721 0.354167361
[24,] 0.6656876 0.668624861 0.334312430
[25,] 0.6926070 0.614786044 0.307393022
[26,] 0.7154596 0.569080730 0.284540365
[27,] 0.7279491 0.544101713 0.272050857
[28,] 0.7369802 0.526039668 0.263019834
[29,] 0.7441034 0.511793222 0.255896611
[30,] 0.7484804 0.503039221 0.251519610
[31,] 0.7508030 0.498393958 0.249196979
[32,] 0.7516135 0.496772993 0.248386497
[33,] 0.7351071 0.529785853 0.264892927
[34,] 0.7396109 0.520778215 0.260389107
[35,] 0.7660249 0.467950274 0.233975137
[36,] 0.7927087 0.414582633 0.207291317
[37,] 0.8575140 0.284971977 0.142485988
[38,] 0.9165802 0.166839571 0.083419786
[39,] 0.9509762 0.098047656 0.049023828
[40,] 0.9758883 0.048223371 0.024111686
[41,] 0.9918608 0.016278376 0.008139188
[42,] 0.9960916 0.007816886 0.003908443
[43,] 0.9978310 0.004337991 0.002168996
[44,] 0.9977603 0.004479464 0.002239732
[45,] 0.9975736 0.004852753 0.002426376
[46,] 0.9964902 0.007019615 0.003509808
[47,] 0.9963543 0.007291417 0.003645709
[48,] 0.9958494 0.008301299 0.004150650
[49,] 0.9970258 0.005948309 0.002974154
[50,] 0.9943242 0.011351549 0.005675775
[51,] 0.9915824 0.016835214 0.008417607
[52,] 0.9908044 0.018391247 0.009195624
[53,] 0.9851396 0.029720894 0.014860447
[54,] 0.9819711 0.036057864 0.018028932
> postscript(file="/var/www/html/rcomp/tmp/1dk2g1258753016.ps",horizontal=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/html/rcomp/tmp/2olhe1258753016.ps",horizontal=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/html/rcomp/tmp/3746x1258753016.ps",horizontal=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/html/rcomp/tmp/4wcko1258753016.ps",horizontal=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/html/rcomp/tmp/5t9zf1258753016.ps",horizontal=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 = 85
Frequency = 1
1 2 3 4 5 6
-70881.1834 -61635.2857 -64098.1429 -61402.7143 -62021.2857 -55116.2857
7 8 9 10 11 12
-56028.4286 -57633.2096 -48765.9239 -43119.2096 -43549.3524 -40035.9239
13 14 15 16 17 18
-28477.1834 -22092.2857 -16222.1429 -15251.7143 -10676.2857 -2640.2857
19 20 21 22 23 24
-1506.4286 -2616.2096 -1706.9239 -5014.2096 -999.3524 3125.0761
25 26 27 28 29 30
15367.8166 23306.7143 18249.8571 20535.2857 25074.7143 21566.7143
31 32 33 34 35 36
29542.5714 32004.7904 43973.0761 41808.7904 39666.6476 31665.0761
37 38 39 40 41 42
43892.8166 52402.7143 47750.8571 51329.2857 51380.7143 48771.7143
43 44 45 46 47 48
47663.5714 47623.7904 44766.0761 43045.7904 41353.6476 34208.0761
49 50 51 52 53 54
44585.8166 47773.7143 48855.8571 56648.2857 55294.7143 47603.7143
55 56 57 58 59 60
44588.5714 39307.7904 20274.0761 13365.7904 1486.6476 1376.0761
61 62 63 64 65 66
3367.8166 -352.2857 569.8571 -5088.7143 -18985.2857 -16601.2857
67 68 69 70 71 72
-16629.4286 -38941.2096 -40280.9239 -42512.2096 -41529.3524 -42277.9239
73 74 75 76 77 78
-37925.1834 -39403.2857 -35106.1429 -46769.7143 -40067.2857 -43584.2857
79 80 81 82 83 84
-47630.4286 -19745.7425 -18259.4567 -7574.7425 3571.1147 11939.5433
85
30069.2837
> postscript(file="/var/www/html/rcomp/tmp/64p6k1258753016.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 85
Frequency = 1
lag(myerror, k = 1) myerror
0 -70881.1834 NA
1 -61635.2857 -70881.1834
2 -64098.1429 -61635.2857
3 -61402.7143 -64098.1429
4 -62021.2857 -61402.7143
5 -55116.2857 -62021.2857
6 -56028.4286 -55116.2857
7 -57633.2096 -56028.4286
8 -48765.9239 -57633.2096
9 -43119.2096 -48765.9239
10 -43549.3524 -43119.2096
11 -40035.9239 -43549.3524
12 -28477.1834 -40035.9239
13 -22092.2857 -28477.1834
14 -16222.1429 -22092.2857
15 -15251.7143 -16222.1429
16 -10676.2857 -15251.7143
17 -2640.2857 -10676.2857
18 -1506.4286 -2640.2857
19 -2616.2096 -1506.4286
20 -1706.9239 -2616.2096
21 -5014.2096 -1706.9239
22 -999.3524 -5014.2096
23 3125.0761 -999.3524
24 15367.8166 3125.0761
25 23306.7143 15367.8166
26 18249.8571 23306.7143
27 20535.2857 18249.8571
28 25074.7143 20535.2857
29 21566.7143 25074.7143
30 29542.5714 21566.7143
31 32004.7904 29542.5714
32 43973.0761 32004.7904
33 41808.7904 43973.0761
34 39666.6476 41808.7904
35 31665.0761 39666.6476
36 43892.8166 31665.0761
37 52402.7143 43892.8166
38 47750.8571 52402.7143
39 51329.2857 47750.8571
40 51380.7143 51329.2857
41 48771.7143 51380.7143
42 47663.5714 48771.7143
43 47623.7904 47663.5714
44 44766.0761 47623.7904
45 43045.7904 44766.0761
46 41353.6476 43045.7904
47 34208.0761 41353.6476
48 44585.8166 34208.0761
49 47773.7143 44585.8166
50 48855.8571 47773.7143
51 56648.2857 48855.8571
52 55294.7143 56648.2857
53 47603.7143 55294.7143
54 44588.5714 47603.7143
55 39307.7904 44588.5714
56 20274.0761 39307.7904
57 13365.7904 20274.0761
58 1486.6476 13365.7904
59 1376.0761 1486.6476
60 3367.8166 1376.0761
61 -352.2857 3367.8166
62 569.8571 -352.2857
63 -5088.7143 569.8571
64 -18985.2857 -5088.7143
65 -16601.2857 -18985.2857
66 -16629.4286 -16601.2857
67 -38941.2096 -16629.4286
68 -40280.9239 -38941.2096
69 -42512.2096 -40280.9239
70 -41529.3524 -42512.2096
71 -42277.9239 -41529.3524
72 -37925.1834 -42277.9239
73 -39403.2857 -37925.1834
74 -35106.1429 -39403.2857
75 -46769.7143 -35106.1429
76 -40067.2857 -46769.7143
77 -43584.2857 -40067.2857
78 -47630.4286 -43584.2857
79 -19745.7425 -47630.4286
80 -18259.4567 -19745.7425
81 -7574.7425 -18259.4567
82 3571.1147 -7574.7425
83 11939.5433 3571.1147
84 30069.2837 11939.5433
85 NA 30069.2837
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -61635.2857 -70881.1834
[2,] -64098.1429 -61635.2857
[3,] -61402.7143 -64098.1429
[4,] -62021.2857 -61402.7143
[5,] -55116.2857 -62021.2857
[6,] -56028.4286 -55116.2857
[7,] -57633.2096 -56028.4286
[8,] -48765.9239 -57633.2096
[9,] -43119.2096 -48765.9239
[10,] -43549.3524 -43119.2096
[11,] -40035.9239 -43549.3524
[12,] -28477.1834 -40035.9239
[13,] -22092.2857 -28477.1834
[14,] -16222.1429 -22092.2857
[15,] -15251.7143 -16222.1429
[16,] -10676.2857 -15251.7143
[17,] -2640.2857 -10676.2857
[18,] -1506.4286 -2640.2857
[19,] -2616.2096 -1506.4286
[20,] -1706.9239 -2616.2096
[21,] -5014.2096 -1706.9239
[22,] -999.3524 -5014.2096
[23,] 3125.0761 -999.3524
[24,] 15367.8166 3125.0761
[25,] 23306.7143 15367.8166
[26,] 18249.8571 23306.7143
[27,] 20535.2857 18249.8571
[28,] 25074.7143 20535.2857
[29,] 21566.7143 25074.7143
[30,] 29542.5714 21566.7143
[31,] 32004.7904 29542.5714
[32,] 43973.0761 32004.7904
[33,] 41808.7904 43973.0761
[34,] 39666.6476 41808.7904
[35,] 31665.0761 39666.6476
[36,] 43892.8166 31665.0761
[37,] 52402.7143 43892.8166
[38,] 47750.8571 52402.7143
[39,] 51329.2857 47750.8571
[40,] 51380.7143 51329.2857
[41,] 48771.7143 51380.7143
[42,] 47663.5714 48771.7143
[43,] 47623.7904 47663.5714
[44,] 44766.0761 47623.7904
[45,] 43045.7904 44766.0761
[46,] 41353.6476 43045.7904
[47,] 34208.0761 41353.6476
[48,] 44585.8166 34208.0761
[49,] 47773.7143 44585.8166
[50,] 48855.8571 47773.7143
[51,] 56648.2857 48855.8571
[52,] 55294.7143 56648.2857
[53,] 47603.7143 55294.7143
[54,] 44588.5714 47603.7143
[55,] 39307.7904 44588.5714
[56,] 20274.0761 39307.7904
[57,] 13365.7904 20274.0761
[58,] 1486.6476 13365.7904
[59,] 1376.0761 1486.6476
[60,] 3367.8166 1376.0761
[61,] -352.2857 3367.8166
[62,] 569.8571 -352.2857
[63,] -5088.7143 569.8571
[64,] -18985.2857 -5088.7143
[65,] -16601.2857 -18985.2857
[66,] -16629.4286 -16601.2857
[67,] -38941.2096 -16629.4286
[68,] -40280.9239 -38941.2096
[69,] -42512.2096 -40280.9239
[70,] -41529.3524 -42512.2096
[71,] -42277.9239 -41529.3524
[72,] -37925.1834 -42277.9239
[73,] -39403.2857 -37925.1834
[74,] -35106.1429 -39403.2857
[75,] -46769.7143 -35106.1429
[76,] -40067.2857 -46769.7143
[77,] -43584.2857 -40067.2857
[78,] -47630.4286 -43584.2857
[79,] -19745.7425 -47630.4286
[80,] -18259.4567 -19745.7425
[81,] -7574.7425 -18259.4567
[82,] 3571.1147 -7574.7425
[83,] 11939.5433 3571.1147
[84,] 30069.2837 11939.5433
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -61635.2857 -70881.1834
2 -64098.1429 -61635.2857
3 -61402.7143 -64098.1429
4 -62021.2857 -61402.7143
5 -55116.2857 -62021.2857
6 -56028.4286 -55116.2857
7 -57633.2096 -56028.4286
8 -48765.9239 -57633.2096
9 -43119.2096 -48765.9239
10 -43549.3524 -43119.2096
11 -40035.9239 -43549.3524
12 -28477.1834 -40035.9239
13 -22092.2857 -28477.1834
14 -16222.1429 -22092.2857
15 -15251.7143 -16222.1429
16 -10676.2857 -15251.7143
17 -2640.2857 -10676.2857
18 -1506.4286 -2640.2857
19 -2616.2096 -1506.4286
20 -1706.9239 -2616.2096
21 -5014.2096 -1706.9239
22 -999.3524 -5014.2096
23 3125.0761 -999.3524
24 15367.8166 3125.0761
25 23306.7143 15367.8166
26 18249.8571 23306.7143
27 20535.2857 18249.8571
28 25074.7143 20535.2857
29 21566.7143 25074.7143
30 29542.5714 21566.7143
31 32004.7904 29542.5714
32 43973.0761 32004.7904
33 41808.7904 43973.0761
34 39666.6476 41808.7904
35 31665.0761 39666.6476
36 43892.8166 31665.0761
37 52402.7143 43892.8166
38 47750.8571 52402.7143
39 51329.2857 47750.8571
40 51380.7143 51329.2857
41 48771.7143 51380.7143
42 47663.5714 48771.7143
43 47623.7904 47663.5714
44 44766.0761 47623.7904
45 43045.7904 44766.0761
46 41353.6476 43045.7904
47 34208.0761 41353.6476
48 44585.8166 34208.0761
49 47773.7143 44585.8166
50 48855.8571 47773.7143
51 56648.2857 48855.8571
52 55294.7143 56648.2857
53 47603.7143 55294.7143
54 44588.5714 47603.7143
55 39307.7904 44588.5714
56 20274.0761 39307.7904
57 13365.7904 20274.0761
58 1486.6476 13365.7904
59 1376.0761 1486.6476
60 3367.8166 1376.0761
61 -352.2857 3367.8166
62 569.8571 -352.2857
63 -5088.7143 569.8571
64 -18985.2857 -5088.7143
65 -16601.2857 -18985.2857
66 -16629.4286 -16601.2857
67 -38941.2096 -16629.4286
68 -40280.9239 -38941.2096
69 -42512.2096 -40280.9239
70 -41529.3524 -42512.2096
71 -42277.9239 -41529.3524
72 -37925.1834 -42277.9239
73 -39403.2857 -37925.1834
74 -35106.1429 -39403.2857
75 -46769.7143 -35106.1429
76 -40067.2857 -46769.7143
77 -43584.2857 -40067.2857
78 -47630.4286 -43584.2857
79 -19745.7425 -47630.4286
80 -18259.4567 -19745.7425
81 -7574.7425 -18259.4567
82 3571.1147 -7574.7425
83 11939.5433 3571.1147
84 30069.2837 11939.5433
> 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/html/rcomp/tmp/7broa1258753016.ps",horizontal=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/html/rcomp/tmp/8b5nm1258753016.ps",horizontal=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/html/rcomp/tmp/9mc3y1258753016.ps",horizontal=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/html/rcomp/tmp/102oqq1258753016.ps",horizontal=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/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/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/html/rcomp/tmp/114xmy1258753016.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/html/rcomp/tmp/12m8xe1258753016.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/html/rcomp/tmp/13shcy1258753017.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/html/rcomp/tmp/14fi8c1258753017.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/html/rcomp/tmp/15nwr31258753017.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/html/rcomp/tmp/16gaqn1258753017.tab")
+ }
>
> system("convert tmp/1dk2g1258753016.ps tmp/1dk2g1258753016.png")
> system("convert tmp/2olhe1258753016.ps tmp/2olhe1258753016.png")
> system("convert tmp/3746x1258753016.ps tmp/3746x1258753016.png")
> system("convert tmp/4wcko1258753016.ps tmp/4wcko1258753016.png")
> system("convert tmp/5t9zf1258753016.ps tmp/5t9zf1258753016.png")
> system("convert tmp/64p6k1258753016.ps tmp/64p6k1258753016.png")
> system("convert tmp/7broa1258753016.ps tmp/7broa1258753016.png")
> system("convert tmp/8b5nm1258753016.ps tmp/8b5nm1258753016.png")
> system("convert tmp/9mc3y1258753016.ps tmp/9mc3y1258753016.png")
> system("convert tmp/102oqq1258753016.ps tmp/102oqq1258753016.png")
>
>
> proc.time()
user system elapsed
2.698 1.563 3.172