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(100
+ ,100
+ ,97.82226485
+ ,99.87129987
+ ,94.04971502
+ ,99.54459954
+ ,91.12460521
+ ,99.81189981
+ ,93.13202153
+ ,100.4851005
+ ,93.88342812
+ ,101.1385011
+ ,92.55349954
+ ,101.3662014
+ ,94.43494835
+ ,101.5147015
+ ,96.25017563
+ ,101.8216018
+ ,100.4355715
+ ,102.4354024
+ ,101.5036685
+ ,102.5344025
+ ,99.39789728
+ ,102.6532027
+ ,99.68990733
+ ,102.4651025
+ ,101.6895041
+ ,102.4354024
+ ,103.6652759
+ ,102.4156024
+ ,103.0532766
+ ,102.4453024
+ ,100.9500712
+ ,102.8908029
+ ,102.345366
+ ,102.8512029
+ ,101.6472299
+ ,103.3561034
+ ,99.56809393
+ ,103.7422037
+ ,95.67727392
+ ,103.7224037
+ ,96.58494865
+ ,104.0788041
+ ,96.32604937
+ ,104.2075042
+ ,95.37109101
+ ,103.9105039
+ ,96.00056203
+ ,103.7026037
+ ,96.88367859
+ ,103.960004
+ ,94.85280372
+ ,104.0986041
+ ,92.46943974
+ ,104.1481041
+ ,93.99180173
+ ,104.7124047
+ ,93.45262168
+ ,104.7223047
+ ,92.26698759
+ ,105.1975052
+ ,90.39653498
+ ,105.0688051
+ ,90.43001228
+ ,105.0589051
+ ,91.04995327
+ ,105.5044055
+ ,89.07845784
+ ,105.3757054
+ ,89.69314509
+ ,105.4747055
+ ,87.92459054
+ ,106.029106
+ ,85.8789319
+ ,107.019107
+ ,83.20612366
+ ,107.3161073
+ ,83.85722053
+ ,107.7517078
+ ,83.01393462
+ ,108.5239085
+ ,82.84508195
+ ,109.3159093
+ ,78.68864276
+ ,109.5634096
+ ,77.56959675
+ ,110.5435105
+ ,78.53689529
+ ,111.1573112
+ ,78.55717715
+ ,111.7414117
+ ,77.4761291
+ ,111.0583111
+ ,81.58931659
+ ,111.2365112
+ ,85.02428326
+ ,111.038511
+ ,91.71290159
+ ,110.3752104
+ ,95.96293061
+ ,110.1376101
+ ,90.84689022
+ ,110.2465102
+ ,92.28788036
+ ,110.6227106
+ ,95.56511274
+ ,109.98911
+ ,93.62452884
+ ,110.2168102
+ ,92.63071726
+ ,110.1376101
+ ,89.50914211
+ ,109.9297099
+ ,87.17171779
+ ,109.8604099
+ ,86.72624975
+ ,110.1970102
+ ,85.63212844
+ ,109.9099099)
+ ,dim=c(2
+ ,60)
+ ,dimnames=list(c('wisselkoers'
+ ,'consumptieprijzen')
+ ,1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('wisselkoers','consumptieprijzen'),1:60))
> 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
wisselkoers consumptieprijzen M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 100.00000 100.0000 1 0 0 0 0 0 0 0 0 0 0
2 97.82226 99.8713 0 1 0 0 0 0 0 0 0 0 0
3 94.04972 99.5446 0 0 1 0 0 0 0 0 0 0 0
4 91.12461 99.8119 0 0 0 1 0 0 0 0 0 0 0
5 93.13202 100.4851 0 0 0 0 1 0 0 0 0 0 0
6 93.88343 101.1385 0 0 0 0 0 1 0 0 0 0 0
7 92.55350 101.3662 0 0 0 0 0 0 1 0 0 0 0
8 94.43495 101.5147 0 0 0 0 0 0 0 1 0 0 0
9 96.25018 101.8216 0 0 0 0 0 0 0 0 1 0 0
10 100.43557 102.4354 0 0 0 0 0 0 0 0 0 1 0
11 101.50367 102.5344 0 0 0 0 0 0 0 0 0 0 1
12 99.39790 102.6532 0 0 0 0 0 0 0 0 0 0 0
13 99.68991 102.4651 1 0 0 0 0 0 0 0 0 0 0
14 101.68950 102.4354 0 1 0 0 0 0 0 0 0 0 0
15 103.66528 102.4156 0 0 1 0 0 0 0 0 0 0 0
16 103.05328 102.4453 0 0 0 1 0 0 0 0 0 0 0
17 100.95007 102.8908 0 0 0 0 1 0 0 0 0 0 0
18 102.34537 102.8512 0 0 0 0 0 1 0 0 0 0 0
19 101.64723 103.3561 0 0 0 0 0 0 1 0 0 0 0
20 99.56809 103.7422 0 0 0 0 0 0 0 1 0 0 0
21 95.67727 103.7224 0 0 0 0 0 0 0 0 1 0 0
22 96.58495 104.0788 0 0 0 0 0 0 0 0 0 1 0
23 96.32605 104.2075 0 0 0 0 0 0 0 0 0 0 1
24 95.37109 103.9105 0 0 0 0 0 0 0 0 0 0 0
25 96.00056 103.7026 1 0 0 0 0 0 0 0 0 0 0
26 96.88368 103.9600 0 1 0 0 0 0 0 0 0 0 0
27 94.85280 104.0986 0 0 1 0 0 0 0 0 0 0 0
28 92.46944 104.1481 0 0 0 1 0 0 0 0 0 0 0
29 93.99180 104.7124 0 0 0 0 1 0 0 0 0 0 0
30 93.45262 104.7223 0 0 0 0 0 1 0 0 0 0 0
31 92.26699 105.1975 0 0 0 0 0 0 1 0 0 0 0
32 90.39653 105.0688 0 0 0 0 0 0 0 1 0 0 0
33 90.43001 105.0589 0 0 0 0 0 0 0 0 1 0 0
34 91.04995 105.5044 0 0 0 0 0 0 0 0 0 1 0
35 89.07846 105.3757 0 0 0 0 0 0 0 0 0 0 1
36 89.69315 105.4747 0 0 0 0 0 0 0 0 0 0 0
37 87.92459 106.0291 1 0 0 0 0 0 0 0 0 0 0
38 85.87893 107.0191 0 1 0 0 0 0 0 0 0 0 0
39 83.20612 107.3161 0 0 1 0 0 0 0 0 0 0 0
40 83.85722 107.7517 0 0 0 1 0 0 0 0 0 0 0
41 83.01393 108.5239 0 0 0 0 1 0 0 0 0 0 0
42 82.84508 109.3159 0 0 0 0 0 1 0 0 0 0 0
43 78.68864 109.5634 0 0 0 0 0 0 1 0 0 0 0
44 77.56960 110.5435 0 0 0 0 0 0 0 1 0 0 0
45 78.53690 111.1573 0 0 0 0 0 0 0 0 1 0 0
46 78.55718 111.7414 0 0 0 0 0 0 0 0 0 1 0
47 77.47613 111.0583 0 0 0 0 0 0 0 0 0 0 1
48 81.58932 111.2365 0 0 0 0 0 0 0 0 0 0 0
49 85.02428 111.0385 1 0 0 0 0 0 0 0 0 0 0
50 91.71290 110.3752 0 1 0 0 0 0 0 0 0 0 0
51 95.96293 110.1376 0 0 1 0 0 0 0 0 0 0 0
52 90.84689 110.2465 0 0 0 1 0 0 0 0 0 0 0
53 92.28788 110.6227 0 0 0 0 1 0 0 0 0 0 0
54 95.56511 109.9891 0 0 0 0 0 1 0 0 0 0 0
55 93.62453 110.2168 0 0 0 0 0 0 1 0 0 0 0
56 92.63072 110.1376 0 0 0 0 0 0 0 1 0 0 0
57 89.50914 109.9297 0 0 0 0 0 0 0 0 1 0 0
58 87.17172 109.8604 0 0 0 0 0 0 0 0 0 1 0
59 86.72625 110.1970 0 0 0 0 0 0 0 0 0 0 1
60 85.63213 109.9099 0 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) consumptieprijzen M1 M2
229.99678 -1.30968 0.78502 1.96612
M3 M4 M5 M6
1.47713 -0.36657 0.77993 1.92797
M7 M8 M9 M10
0.50667 0.01277 -0.64758 0.53726
M11
-0.06533
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-8.32204 -4.15818 -0.04975 4.07729 8.73380
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 229.99678 20.92982 10.989 1.40e-14 ***
consumptieprijzen -1.30968 0.19496 -6.718 2.19e-08 ***
M1 0.78502 3.43633 0.228 0.820
M2 1.96612 3.43450 0.572 0.570
M3 1.47713 3.43513 0.430 0.669
M4 -0.36657 3.43149 -0.107 0.915
M5 0.77993 3.42223 0.228 0.821
M6 1.92797 3.42030 0.564 0.576
M7 0.50667 3.41706 0.148 0.883
M8 0.01277 3.41542 0.004 0.997
M9 -0.64758 3.41486 -0.190 0.850
M10 0.53726 3.41440 0.157 0.876
M11 -0.06533 3.41437 -0.019 0.985
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.399 on 47 degrees of freedom
Multiple R-squared: 0.5182, Adjusted R-squared: 0.3952
F-statistic: 4.213 on 12 and 47 DF, p-value: 0.0001737
> 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.29925608 0.59851217 0.7007439
[2,] 0.16693402 0.33386804 0.8330660
[3,] 0.11210985 0.22421971 0.8878901
[4,] 0.07424823 0.14849646 0.9257518
[5,] 0.04006775 0.08013549 0.9599323
[6,] 0.04032752 0.08065503 0.9596725
[7,] 0.06344235 0.12688470 0.9365577
[8,] 0.10362812 0.20725624 0.8963719
[9,] 0.09983238 0.19966476 0.9001676
[10,] 0.13237010 0.26474019 0.8676299
[11,] 0.13037067 0.26074135 0.8696293
[12,] 0.12638394 0.25276789 0.8736161
[13,] 0.11364542 0.22729084 0.8863546
[14,] 0.08488327 0.16976654 0.9151167
[15,] 0.06567490 0.13134981 0.9343251
[16,] 0.04974786 0.09949572 0.9502521
[17,] 0.04087474 0.08174947 0.9591253
[18,] 0.02973058 0.05946117 0.9702694
[19,] 0.02734333 0.05468666 0.9726567
[20,] 0.03159185 0.06318371 0.9684081
[21,] 0.03104418 0.06208837 0.9689558
[22,] 0.03307948 0.06615895 0.9669205
[23,] 0.02604959 0.05209919 0.9739504
[24,] 0.02863977 0.05727955 0.9713602
[25,] 0.02040136 0.04080273 0.9795986
[26,] 0.02815048 0.05630097 0.9718495
[27,] 0.05976386 0.11952771 0.9402361
[28,] 0.35666481 0.71332962 0.6433352
[29,] 0.77006924 0.45986152 0.2299308
> postscript(file="/var/www/html/rcomp/tmp/1fq9o1258734167.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/2eug21258734167.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/323te1258734167.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/4otp91258734167.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/52v2m1258734167.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 = 60
Frequency = 1
1 2 3 4 5 6 7
0.1859742 -3.3414103 -7.0528484 -7.7841824 -6.0415907 -5.5824800 -5.1928882
8 9 10 11 12 13 14
-2.6230503 0.2544675 4.0589034 5.8592505 3.8437346 3.1043714 3.8839770
15 16 17 18 19 20 21
6.3228011 7.5933977 4.9271539 5.1225453 6.5069725 5.4274053 2.1710037
22 23 24 25 26 27 28
2.3606072 2.8728554 1.4635877 1.0357539 1.0748883 -0.2854812 -0.7603176
29 30 31 32 33 34 35
0.3545958 -1.3196587 -0.4616269 -2.0067333 -1.3258718 -1.3073098 -2.8447690
36 37 38 39 40 41 42
-2.1657581 -3.9932493 -5.9234193 -7.7182689 -4.6529773 -5.6314296 -5.9110567
43 44 45 46 47 48 49
-8.3220438 -7.6635717 -5.2320420 -5.6316176 -7.0047155 -2.7234779 -0.3328503
50 51 52 53 54 55 56
4.3059643 8.7337974 5.6040796 6.3912705 7.6906501 7.4695865 6.8659501
57 58 59 60
4.1324427 0.5194168 1.1173785 -0.4180863
> postscript(file="/var/www/html/rcomp/tmp/69y2w1258734167.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 0.1859742 NA
1 -3.3414103 0.1859742
2 -7.0528484 -3.3414103
3 -7.7841824 -7.0528484
4 -6.0415907 -7.7841824
5 -5.5824800 -6.0415907
6 -5.1928882 -5.5824800
7 -2.6230503 -5.1928882
8 0.2544675 -2.6230503
9 4.0589034 0.2544675
10 5.8592505 4.0589034
11 3.8437346 5.8592505
12 3.1043714 3.8437346
13 3.8839770 3.1043714
14 6.3228011 3.8839770
15 7.5933977 6.3228011
16 4.9271539 7.5933977
17 5.1225453 4.9271539
18 6.5069725 5.1225453
19 5.4274053 6.5069725
20 2.1710037 5.4274053
21 2.3606072 2.1710037
22 2.8728554 2.3606072
23 1.4635877 2.8728554
24 1.0357539 1.4635877
25 1.0748883 1.0357539
26 -0.2854812 1.0748883
27 -0.7603176 -0.2854812
28 0.3545958 -0.7603176
29 -1.3196587 0.3545958
30 -0.4616269 -1.3196587
31 -2.0067333 -0.4616269
32 -1.3258718 -2.0067333
33 -1.3073098 -1.3258718
34 -2.8447690 -1.3073098
35 -2.1657581 -2.8447690
36 -3.9932493 -2.1657581
37 -5.9234193 -3.9932493
38 -7.7182689 -5.9234193
39 -4.6529773 -7.7182689
40 -5.6314296 -4.6529773
41 -5.9110567 -5.6314296
42 -8.3220438 -5.9110567
43 -7.6635717 -8.3220438
44 -5.2320420 -7.6635717
45 -5.6316176 -5.2320420
46 -7.0047155 -5.6316176
47 -2.7234779 -7.0047155
48 -0.3328503 -2.7234779
49 4.3059643 -0.3328503
50 8.7337974 4.3059643
51 5.6040796 8.7337974
52 6.3912705 5.6040796
53 7.6906501 6.3912705
54 7.4695865 7.6906501
55 6.8659501 7.4695865
56 4.1324427 6.8659501
57 0.5194168 4.1324427
58 1.1173785 0.5194168
59 -0.4180863 1.1173785
60 NA -0.4180863
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -3.3414103 0.1859742
[2,] -7.0528484 -3.3414103
[3,] -7.7841824 -7.0528484
[4,] -6.0415907 -7.7841824
[5,] -5.5824800 -6.0415907
[6,] -5.1928882 -5.5824800
[7,] -2.6230503 -5.1928882
[8,] 0.2544675 -2.6230503
[9,] 4.0589034 0.2544675
[10,] 5.8592505 4.0589034
[11,] 3.8437346 5.8592505
[12,] 3.1043714 3.8437346
[13,] 3.8839770 3.1043714
[14,] 6.3228011 3.8839770
[15,] 7.5933977 6.3228011
[16,] 4.9271539 7.5933977
[17,] 5.1225453 4.9271539
[18,] 6.5069725 5.1225453
[19,] 5.4274053 6.5069725
[20,] 2.1710037 5.4274053
[21,] 2.3606072 2.1710037
[22,] 2.8728554 2.3606072
[23,] 1.4635877 2.8728554
[24,] 1.0357539 1.4635877
[25,] 1.0748883 1.0357539
[26,] -0.2854812 1.0748883
[27,] -0.7603176 -0.2854812
[28,] 0.3545958 -0.7603176
[29,] -1.3196587 0.3545958
[30,] -0.4616269 -1.3196587
[31,] -2.0067333 -0.4616269
[32,] -1.3258718 -2.0067333
[33,] -1.3073098 -1.3258718
[34,] -2.8447690 -1.3073098
[35,] -2.1657581 -2.8447690
[36,] -3.9932493 -2.1657581
[37,] -5.9234193 -3.9932493
[38,] -7.7182689 -5.9234193
[39,] -4.6529773 -7.7182689
[40,] -5.6314296 -4.6529773
[41,] -5.9110567 -5.6314296
[42,] -8.3220438 -5.9110567
[43,] -7.6635717 -8.3220438
[44,] -5.2320420 -7.6635717
[45,] -5.6316176 -5.2320420
[46,] -7.0047155 -5.6316176
[47,] -2.7234779 -7.0047155
[48,] -0.3328503 -2.7234779
[49,] 4.3059643 -0.3328503
[50,] 8.7337974 4.3059643
[51,] 5.6040796 8.7337974
[52,] 6.3912705 5.6040796
[53,] 7.6906501 6.3912705
[54,] 7.4695865 7.6906501
[55,] 6.8659501 7.4695865
[56,] 4.1324427 6.8659501
[57,] 0.5194168 4.1324427
[58,] 1.1173785 0.5194168
[59,] -0.4180863 1.1173785
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -3.3414103 0.1859742
2 -7.0528484 -3.3414103
3 -7.7841824 -7.0528484
4 -6.0415907 -7.7841824
5 -5.5824800 -6.0415907
6 -5.1928882 -5.5824800
7 -2.6230503 -5.1928882
8 0.2544675 -2.6230503
9 4.0589034 0.2544675
10 5.8592505 4.0589034
11 3.8437346 5.8592505
12 3.1043714 3.8437346
13 3.8839770 3.1043714
14 6.3228011 3.8839770
15 7.5933977 6.3228011
16 4.9271539 7.5933977
17 5.1225453 4.9271539
18 6.5069725 5.1225453
19 5.4274053 6.5069725
20 2.1710037 5.4274053
21 2.3606072 2.1710037
22 2.8728554 2.3606072
23 1.4635877 2.8728554
24 1.0357539 1.4635877
25 1.0748883 1.0357539
26 -0.2854812 1.0748883
27 -0.7603176 -0.2854812
28 0.3545958 -0.7603176
29 -1.3196587 0.3545958
30 -0.4616269 -1.3196587
31 -2.0067333 -0.4616269
32 -1.3258718 -2.0067333
33 -1.3073098 -1.3258718
34 -2.8447690 -1.3073098
35 -2.1657581 -2.8447690
36 -3.9932493 -2.1657581
37 -5.9234193 -3.9932493
38 -7.7182689 -5.9234193
39 -4.6529773 -7.7182689
40 -5.6314296 -4.6529773
41 -5.9110567 -5.6314296
42 -8.3220438 -5.9110567
43 -7.6635717 -8.3220438
44 -5.2320420 -7.6635717
45 -5.6316176 -5.2320420
46 -7.0047155 -5.6316176
47 -2.7234779 -7.0047155
48 -0.3328503 -2.7234779
49 4.3059643 -0.3328503
50 8.7337974 4.3059643
51 5.6040796 8.7337974
52 6.3912705 5.6040796
53 7.6906501 6.3912705
54 7.4695865 7.6906501
55 6.8659501 7.4695865
56 4.1324427 6.8659501
57 0.5194168 4.1324427
58 1.1173785 0.5194168
59 -0.4180863 1.1173785
> 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/7ecye1258734167.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/854if1258734167.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/9kvh61258734167.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/108vdv1258734167.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/11jgz31258734167.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/12ufve1258734167.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/13vixn1258734167.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/14yaaw1258734167.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/15xfwe1258734167.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/16zda91258734168.tab")
+ }
>
> system("convert tmp/1fq9o1258734167.ps tmp/1fq9o1258734167.png")
> system("convert tmp/2eug21258734167.ps tmp/2eug21258734167.png")
> system("convert tmp/323te1258734167.ps tmp/323te1258734167.png")
> system("convert tmp/4otp91258734167.ps tmp/4otp91258734167.png")
> system("convert tmp/52v2m1258734167.ps tmp/52v2m1258734167.png")
> system("convert tmp/69y2w1258734167.ps tmp/69y2w1258734167.png")
> system("convert tmp/7ecye1258734167.ps tmp/7ecye1258734167.png")
> system("convert tmp/854if1258734167.ps tmp/854if1258734167.png")
> system("convert tmp/9kvh61258734167.ps tmp/9kvh61258734167.png")
> system("convert tmp/108vdv1258734167.ps tmp/108vdv1258734167.png")
>
>
> proc.time()
user system elapsed
2.397 1.564 2.771