R version 2.8.0 (2008-10-20)
Copyright (C) 2008 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(104.3,0,103.9,0,103.9,0,103.9,0,108.0,0,108.0,0,108.0,0,108.0,0,108.0,0,108.0,0,108.0,0,108.0,0,108.0,0,108.0,0,108.0,0,108.0,0,108.0,0,108.0,0,108.0,0,108.0,0,108.0,0,108.0,0,108.0,1,108.0,1,108.0,1,108.0,1,108.0,1,108.2,1,112.3,1,111.3,1,111.3,1,115.3,1,117.2,1,118.3,1,118.3,1,118.3,1,119.0,1,120.6,1,122.6,1,122.6,1,127.4,1,125.9,1,121.5,1,118.8,1,121.6,1,122.3,1,122.7,1,120.8,1,120.1,1,120.1,1,120.1,1,120.1,1,128.4,1,129.8,1,129.8,1,128.6,1,128.6,1,133.7,1,130.0,1,125.9,1,129.4,1,129.4,1,130.6,1,130.6,1,130.6,1,130.8,1,129.7,1,125.8,1,126.0,1,125.6,1,125.4,1,124.7,1,126.9,1,129.1,1),dim=c(2,74),dimnames=list(c('y','x'),1:74))
> y <- array(NA,dim=c(2,74),dimnames=list(c('y','x'),1:74))
> 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'
> #'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
y x M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 104.3 0 1 0 0 0 0 0 0 0 0 0 0 1
2 103.9 0 0 1 0 0 0 0 0 0 0 0 0 2
3 103.9 0 0 0 1 0 0 0 0 0 0 0 0 3
4 103.9 0 0 0 0 1 0 0 0 0 0 0 0 4
5 108.0 0 0 0 0 0 1 0 0 0 0 0 0 5
6 108.0 0 0 0 0 0 0 1 0 0 0 0 0 6
7 108.0 0 0 0 0 0 0 0 1 0 0 0 0 7
8 108.0 0 0 0 0 0 0 0 0 1 0 0 0 8
9 108.0 0 0 0 0 0 0 0 0 0 1 0 0 9
10 108.0 0 0 0 0 0 0 0 0 0 0 1 0 10
11 108.0 0 0 0 0 0 0 0 0 0 0 0 1 11
12 108.0 0 0 0 0 0 0 0 0 0 0 0 0 12
13 108.0 0 1 0 0 0 0 0 0 0 0 0 0 13
14 108.0 0 0 1 0 0 0 0 0 0 0 0 0 14
15 108.0 0 0 0 1 0 0 0 0 0 0 0 0 15
16 108.0 0 0 0 0 1 0 0 0 0 0 0 0 16
17 108.0 0 0 0 0 0 1 0 0 0 0 0 0 17
18 108.0 0 0 0 0 0 0 1 0 0 0 0 0 18
19 108.0 0 0 0 0 0 0 0 1 0 0 0 0 19
20 108.0 0 0 0 0 0 0 0 0 1 0 0 0 20
21 108.0 0 0 0 0 0 0 0 0 0 1 0 0 21
22 108.0 0 0 0 0 0 0 0 0 0 0 1 0 22
23 108.0 1 0 0 0 0 0 0 0 0 0 0 1 23
24 108.0 1 0 0 0 0 0 0 0 0 0 0 0 24
25 108.0 1 1 0 0 0 0 0 0 0 0 0 0 25
26 108.0 1 0 1 0 0 0 0 0 0 0 0 0 26
27 108.0 1 0 0 1 0 0 0 0 0 0 0 0 27
28 108.2 1 0 0 0 1 0 0 0 0 0 0 0 28
29 112.3 1 0 0 0 0 1 0 0 0 0 0 0 29
30 111.3 1 0 0 0 0 0 1 0 0 0 0 0 30
31 111.3 1 0 0 0 0 0 0 1 0 0 0 0 31
32 115.3 1 0 0 0 0 0 0 0 1 0 0 0 32
33 117.2 1 0 0 0 0 0 0 0 0 1 0 0 33
34 118.3 1 0 0 0 0 0 0 0 0 0 1 0 34
35 118.3 1 0 0 0 0 0 0 0 0 0 0 1 35
36 118.3 1 0 0 0 0 0 0 0 0 0 0 0 36
37 119.0 1 1 0 0 0 0 0 0 0 0 0 0 37
38 120.6 1 0 1 0 0 0 0 0 0 0 0 0 38
39 122.6 1 0 0 1 0 0 0 0 0 0 0 0 39
40 122.6 1 0 0 0 1 0 0 0 0 0 0 0 40
41 127.4 1 0 0 0 0 1 0 0 0 0 0 0 41
42 125.9 1 0 0 0 0 0 1 0 0 0 0 0 42
43 121.5 1 0 0 0 0 0 0 1 0 0 0 0 43
44 118.8 1 0 0 0 0 0 0 0 1 0 0 0 44
45 121.6 1 0 0 0 0 0 0 0 0 1 0 0 45
46 122.3 1 0 0 0 0 0 0 0 0 0 1 0 46
47 122.7 1 0 0 0 0 0 0 0 0 0 0 1 47
48 120.8 1 0 0 0 0 0 0 0 0 0 0 0 48
49 120.1 1 1 0 0 0 0 0 0 0 0 0 0 49
50 120.1 1 0 1 0 0 0 0 0 0 0 0 0 50
51 120.1 1 0 0 1 0 0 0 0 0 0 0 0 51
52 120.1 1 0 0 0 1 0 0 0 0 0 0 0 52
53 128.4 1 0 0 0 0 1 0 0 0 0 0 0 53
54 129.8 1 0 0 0 0 0 1 0 0 0 0 0 54
55 129.8 1 0 0 0 0 0 0 1 0 0 0 0 55
56 128.6 1 0 0 0 0 0 0 0 1 0 0 0 56
57 128.6 1 0 0 0 0 0 0 0 0 1 0 0 57
58 133.7 1 0 0 0 0 0 0 0 0 0 1 0 58
59 130.0 1 0 0 0 0 0 0 0 0 0 0 1 59
60 125.9 1 0 0 0 0 0 0 0 0 0 0 0 60
61 129.4 1 1 0 0 0 0 0 0 0 0 0 0 61
62 129.4 1 0 1 0 0 0 0 0 0 0 0 0 62
63 130.6 1 0 0 1 0 0 0 0 0 0 0 0 63
64 130.6 1 0 0 0 1 0 0 0 0 0 0 0 64
65 130.6 1 0 0 0 0 1 0 0 0 0 0 0 65
66 130.8 1 0 0 0 0 0 1 0 0 0 0 0 66
67 129.7 1 0 0 0 0 0 0 1 0 0 0 0 67
68 125.8 1 0 0 0 0 0 0 0 1 0 0 0 68
69 126.0 1 0 0 0 0 0 0 0 0 1 0 0 69
70 125.6 1 0 0 0 0 0 0 0 0 0 1 0 70
71 125.4 1 0 0 0 0 0 0 0 0 0 0 1 71
72 124.7 1 0 0 0 0 0 0 0 0 0 0 0 72
73 126.9 1 1 0 0 0 0 0 0 0 0 0 0 73
74 129.1 1 0 1 0 0 0 0 0 0 0 0 0 74
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) x M1 M2 M3 M4
100.7604 0.3702 0.9259 1.0177 1.5243 1.1636
M5 M6 M7 M8 M9 M10
4.3197 3.7757 2.4650 1.4377 1.8603 2.5497
M11 t
1.5107 0.3940
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-5.660 -2.704 0.728 2.329 7.168
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 100.76044 1.65031 61.056 <2e-16 ***
x 0.37015 1.47544 0.251 0.8028
M1 0.92594 1.95032 0.475 0.6367
M2 1.01766 1.95028 0.522 0.6037
M3 1.52431 2.02829 0.752 0.4553
M4 1.16365 2.02715 0.574 0.5681
M5 4.31965 2.02651 2.132 0.0371 *
M6 3.77566 2.02636 1.863 0.0673 .
M7 2.46500 2.02670 1.216 0.2286
M8 1.43767 2.02752 0.709 0.4810
M9 1.86034 2.02884 0.917 0.3628
M10 2.54968 2.03065 1.256 0.2141
M11 1.51066 2.02106 0.747 0.4577
t 0.39399 0.03153 12.496 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.5 on 60 degrees of freedom
Multiple R-squared: 0.8811, Adjusted R-squared: 0.8553
F-statistic: 34.19 on 13 and 60 DF, p-value: < 2.2e-16
> 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,] 6.957980e-02 1.391596e-01 0.9304202
[2,] 5.262039e-02 1.052408e-01 0.9473796
[3,] 3.170326e-02 6.340652e-02 0.9682967
[4,] 1.712558e-02 3.425116e-02 0.9828744
[5,] 8.458695e-03 1.691739e-02 0.9915413
[6,] 3.866482e-03 7.732964e-03 0.9961335
[7,] 1.417546e-03 2.835091e-03 0.9985825
[8,] 4.726891e-04 9.453782e-04 0.9995273
[9,] 1.790066e-04 3.580133e-04 0.9998210
[10,] 7.103063e-05 1.420613e-04 0.9999290
[11,] 3.413323e-05 6.826645e-05 0.9999659
[12,] 1.850568e-05 3.701136e-05 0.9999815
[13,] 3.542740e-05 7.085479e-05 0.9999646
[14,] 5.398684e-05 1.079737e-04 0.9999460
[15,] 7.960516e-05 1.592103e-04 0.9999204
[16,] 7.095692e-04 1.419138e-03 0.9992904
[17,] 6.477514e-03 1.295503e-02 0.9935225
[18,] 2.920258e-02 5.840516e-02 0.9707974
[19,] 6.731196e-02 1.346239e-01 0.9326880
[20,] 9.627796e-02 1.925559e-01 0.9037220
[21,] 1.231964e-01 2.463929e-01 0.8768036
[22,] 1.672569e-01 3.345139e-01 0.8327431
[23,] 2.438657e-01 4.877315e-01 0.7561343
[24,] 2.821489e-01 5.642977e-01 0.7178511
[25,] 3.567501e-01 7.135003e-01 0.6432499
[26,] 3.520765e-01 7.041531e-01 0.6479235
[27,] 3.058959e-01 6.117918e-01 0.6941041
[28,] 2.768932e-01 5.537864e-01 0.7231068
[29,] 2.131711e-01 4.263423e-01 0.7868289
[30,] 1.751467e-01 3.502934e-01 0.8248533
[31,] 1.265315e-01 2.530630e-01 0.8734685
[32,] 9.047703e-02 1.809541e-01 0.9095230
[33,] 9.568685e-02 1.913737e-01 0.9043131
[34,] 1.559615e-01 3.119231e-01 0.8440385
[35,] 3.616957e-01 7.233915e-01 0.6383043
[36,] 8.733168e-01 2.533665e-01 0.1266832
[37,] 8.883879e-01 2.232241e-01 0.1116121
[38,] 8.871773e-01 2.256453e-01 0.1128227
[39,] 8.639345e-01 2.721311e-01 0.1360655
[40,] 7.571509e-01 4.856982e-01 0.2428491
[41,] 5.980360e-01 8.039279e-01 0.4019640
> postscript(file="/var/www/html/freestat/rcomp/tmp/1v01r1230118452.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/freestat/rcomp/tmp/2ujpq1230118452.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/freestat/rcomp/tmp/308sg1230118452.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/freestat/rcomp/tmp/420c21230118452.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/freestat/rcomp/tmp/5g8781230118452.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 = 74
Frequency = 1
1 2 3 4 5 6 7
2.2196264 1.3339121 0.4332711 0.3999377 0.9499377 1.0999377 2.0166044
8 9 10 11 12 13 14
2.6499377 1.8332711 0.7499377 1.3949634 2.5116300 1.1916923 0.7059780
15 16 17 18 19 20 21
-0.1946630 -0.2279963 -3.7779963 -3.6279963 -2.7113297 -2.0779963 -2.8946630
22 23 24 25 26 27 28
-3.9779963 -3.7031245 -2.5864579 -3.9063956 -4.3921099 -5.2927509 -5.1260842
29 30 31 32 33 34 35
-4.5760842 -5.4260842 -4.5094176 0.1239158 1.2072491 1.2239158 1.8689414
36 37 38 39 40 41 42
2.9856081 2.3656703 3.4799560 4.5793150 4.5459817 5.7959817 4.4459817
43 44 45 46 47 48 49
0.9626484 -1.1040183 0.8793150 0.4959817 1.5410073 0.7576740 -1.2622637
50 51 52 53 54 55 56
-1.7479780 -2.6486190 -2.6819524 2.0680476 3.6180476 4.5347143 3.9680476
57 58 59 60 61 62 63
3.1513810 7.1680476 4.1130733 1.1297399 3.3098022 2.8240879 3.1234469
64 65 66 67 68 69 70
3.0901136 -0.4598864 -0.1098864 -0.2932198 -3.5598864 -4.1765531 -5.6598864
71 72 73 74
-5.2148608 -4.7981941 -3.9181319 -2.2038462
> postscript(file="/var/www/html/freestat/rcomp/tmp/64xly1230118452.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 = 74
Frequency = 1
lag(myerror, k = 1) myerror
0 2.2196264 NA
1 1.3339121 2.2196264
2 0.4332711 1.3339121
3 0.3999377 0.4332711
4 0.9499377 0.3999377
5 1.0999377 0.9499377
6 2.0166044 1.0999377
7 2.6499377 2.0166044
8 1.8332711 2.6499377
9 0.7499377 1.8332711
10 1.3949634 0.7499377
11 2.5116300 1.3949634
12 1.1916923 2.5116300
13 0.7059780 1.1916923
14 -0.1946630 0.7059780
15 -0.2279963 -0.1946630
16 -3.7779963 -0.2279963
17 -3.6279963 -3.7779963
18 -2.7113297 -3.6279963
19 -2.0779963 -2.7113297
20 -2.8946630 -2.0779963
21 -3.9779963 -2.8946630
22 -3.7031245 -3.9779963
23 -2.5864579 -3.7031245
24 -3.9063956 -2.5864579
25 -4.3921099 -3.9063956
26 -5.2927509 -4.3921099
27 -5.1260842 -5.2927509
28 -4.5760842 -5.1260842
29 -5.4260842 -4.5760842
30 -4.5094176 -5.4260842
31 0.1239158 -4.5094176
32 1.2072491 0.1239158
33 1.2239158 1.2072491
34 1.8689414 1.2239158
35 2.9856081 1.8689414
36 2.3656703 2.9856081
37 3.4799560 2.3656703
38 4.5793150 3.4799560
39 4.5459817 4.5793150
40 5.7959817 4.5459817
41 4.4459817 5.7959817
42 0.9626484 4.4459817
43 -1.1040183 0.9626484
44 0.8793150 -1.1040183
45 0.4959817 0.8793150
46 1.5410073 0.4959817
47 0.7576740 1.5410073
48 -1.2622637 0.7576740
49 -1.7479780 -1.2622637
50 -2.6486190 -1.7479780
51 -2.6819524 -2.6486190
52 2.0680476 -2.6819524
53 3.6180476 2.0680476
54 4.5347143 3.6180476
55 3.9680476 4.5347143
56 3.1513810 3.9680476
57 7.1680476 3.1513810
58 4.1130733 7.1680476
59 1.1297399 4.1130733
60 3.3098022 1.1297399
61 2.8240879 3.3098022
62 3.1234469 2.8240879
63 3.0901136 3.1234469
64 -0.4598864 3.0901136
65 -0.1098864 -0.4598864
66 -0.2932198 -0.1098864
67 -3.5598864 -0.2932198
68 -4.1765531 -3.5598864
69 -5.6598864 -4.1765531
70 -5.2148608 -5.6598864
71 -4.7981941 -5.2148608
72 -3.9181319 -4.7981941
73 -2.2038462 -3.9181319
74 NA -2.2038462
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 1.3339121 2.2196264
[2,] 0.4332711 1.3339121
[3,] 0.3999377 0.4332711
[4,] 0.9499377 0.3999377
[5,] 1.0999377 0.9499377
[6,] 2.0166044 1.0999377
[7,] 2.6499377 2.0166044
[8,] 1.8332711 2.6499377
[9,] 0.7499377 1.8332711
[10,] 1.3949634 0.7499377
[11,] 2.5116300 1.3949634
[12,] 1.1916923 2.5116300
[13,] 0.7059780 1.1916923
[14,] -0.1946630 0.7059780
[15,] -0.2279963 -0.1946630
[16,] -3.7779963 -0.2279963
[17,] -3.6279963 -3.7779963
[18,] -2.7113297 -3.6279963
[19,] -2.0779963 -2.7113297
[20,] -2.8946630 -2.0779963
[21,] -3.9779963 -2.8946630
[22,] -3.7031245 -3.9779963
[23,] -2.5864579 -3.7031245
[24,] -3.9063956 -2.5864579
[25,] -4.3921099 -3.9063956
[26,] -5.2927509 -4.3921099
[27,] -5.1260842 -5.2927509
[28,] -4.5760842 -5.1260842
[29,] -5.4260842 -4.5760842
[30,] -4.5094176 -5.4260842
[31,] 0.1239158 -4.5094176
[32,] 1.2072491 0.1239158
[33,] 1.2239158 1.2072491
[34,] 1.8689414 1.2239158
[35,] 2.9856081 1.8689414
[36,] 2.3656703 2.9856081
[37,] 3.4799560 2.3656703
[38,] 4.5793150 3.4799560
[39,] 4.5459817 4.5793150
[40,] 5.7959817 4.5459817
[41,] 4.4459817 5.7959817
[42,] 0.9626484 4.4459817
[43,] -1.1040183 0.9626484
[44,] 0.8793150 -1.1040183
[45,] 0.4959817 0.8793150
[46,] 1.5410073 0.4959817
[47,] 0.7576740 1.5410073
[48,] -1.2622637 0.7576740
[49,] -1.7479780 -1.2622637
[50,] -2.6486190 -1.7479780
[51,] -2.6819524 -2.6486190
[52,] 2.0680476 -2.6819524
[53,] 3.6180476 2.0680476
[54,] 4.5347143 3.6180476
[55,] 3.9680476 4.5347143
[56,] 3.1513810 3.9680476
[57,] 7.1680476 3.1513810
[58,] 4.1130733 7.1680476
[59,] 1.1297399 4.1130733
[60,] 3.3098022 1.1297399
[61,] 2.8240879 3.3098022
[62,] 3.1234469 2.8240879
[63,] 3.0901136 3.1234469
[64,] -0.4598864 3.0901136
[65,] -0.1098864 -0.4598864
[66,] -0.2932198 -0.1098864
[67,] -3.5598864 -0.2932198
[68,] -4.1765531 -3.5598864
[69,] -5.6598864 -4.1765531
[70,] -5.2148608 -5.6598864
[71,] -4.7981941 -5.2148608
[72,] -3.9181319 -4.7981941
[73,] -2.2038462 -3.9181319
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 1.3339121 2.2196264
2 0.4332711 1.3339121
3 0.3999377 0.4332711
4 0.9499377 0.3999377
5 1.0999377 0.9499377
6 2.0166044 1.0999377
7 2.6499377 2.0166044
8 1.8332711 2.6499377
9 0.7499377 1.8332711
10 1.3949634 0.7499377
11 2.5116300 1.3949634
12 1.1916923 2.5116300
13 0.7059780 1.1916923
14 -0.1946630 0.7059780
15 -0.2279963 -0.1946630
16 -3.7779963 -0.2279963
17 -3.6279963 -3.7779963
18 -2.7113297 -3.6279963
19 -2.0779963 -2.7113297
20 -2.8946630 -2.0779963
21 -3.9779963 -2.8946630
22 -3.7031245 -3.9779963
23 -2.5864579 -3.7031245
24 -3.9063956 -2.5864579
25 -4.3921099 -3.9063956
26 -5.2927509 -4.3921099
27 -5.1260842 -5.2927509
28 -4.5760842 -5.1260842
29 -5.4260842 -4.5760842
30 -4.5094176 -5.4260842
31 0.1239158 -4.5094176
32 1.2072491 0.1239158
33 1.2239158 1.2072491
34 1.8689414 1.2239158
35 2.9856081 1.8689414
36 2.3656703 2.9856081
37 3.4799560 2.3656703
38 4.5793150 3.4799560
39 4.5459817 4.5793150
40 5.7959817 4.5459817
41 4.4459817 5.7959817
42 0.9626484 4.4459817
43 -1.1040183 0.9626484
44 0.8793150 -1.1040183
45 0.4959817 0.8793150
46 1.5410073 0.4959817
47 0.7576740 1.5410073
48 -1.2622637 0.7576740
49 -1.7479780 -1.2622637
50 -2.6486190 -1.7479780
51 -2.6819524 -2.6486190
52 2.0680476 -2.6819524
53 3.6180476 2.0680476
54 4.5347143 3.6180476
55 3.9680476 4.5347143
56 3.1513810 3.9680476
57 7.1680476 3.1513810
58 4.1130733 7.1680476
59 1.1297399 4.1130733
60 3.3098022 1.1297399
61 2.8240879 3.3098022
62 3.1234469 2.8240879
63 3.0901136 3.1234469
64 -0.4598864 3.0901136
65 -0.1098864 -0.4598864
66 -0.2932198 -0.1098864
67 -3.5598864 -0.2932198
68 -4.1765531 -3.5598864
69 -5.6598864 -4.1765531
70 -5.2148608 -5.6598864
71 -4.7981941 -5.2148608
72 -3.9181319 -4.7981941
73 -2.2038462 -3.9181319
> 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/freestat/rcomp/tmp/7rfi71230118452.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/freestat/rcomp/tmp/82sfd1230118452.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/freestat/rcomp/tmp/93l6b1230118452.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/freestat/rcomp/tmp/105ghp1230118452.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/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/1172i81230118452.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/freestat/rcomp/tmp/12sax51230118452.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/freestat/rcomp/tmp/1361hj1230118452.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/freestat/rcomp/tmp/14uall1230118452.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/freestat/rcomp/tmp/15anpz1230118452.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/freestat/rcomp/tmp/166q691230118452.tab")
+ }
>
> system("convert tmp/1v01r1230118452.ps tmp/1v01r1230118452.png")
> system("convert tmp/2ujpq1230118452.ps tmp/2ujpq1230118452.png")
> system("convert tmp/308sg1230118452.ps tmp/308sg1230118452.png")
> system("convert tmp/420c21230118452.ps tmp/420c21230118452.png")
> system("convert tmp/5g8781230118452.ps tmp/5g8781230118452.png")
> system("convert tmp/64xly1230118452.ps tmp/64xly1230118452.png")
> system("convert tmp/7rfi71230118452.ps tmp/7rfi71230118452.png")
> system("convert tmp/82sfd1230118452.ps tmp/82sfd1230118452.png")
> system("convert tmp/93l6b1230118452.ps tmp/93l6b1230118452.png")
> system("convert tmp/105ghp1230118452.ps tmp/105ghp1230118452.png")
>
>
> proc.time()
user system elapsed
3.918 2.531 4.815