R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-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(1925
+ ,358
+ ,155
+ ,1580
+ ,375
+ ,172
+ ,1961
+ ,761
+ ,467
+ ,1807
+ ,477
+ ,241
+ ,1526
+ ,547
+ ,294
+ ,1802
+ ,879
+ ,567
+ ,1822
+ ,450
+ ,280
+ ,1125
+ ,462
+ ,225
+ ,1569
+ ,1613
+ ,558
+ ,1829
+ ,854
+ ,342
+ ,1575
+ ,761
+ ,309
+ ,2339
+ ,1521
+ ,1437
+ ,2355
+ ,666
+ ,241
+ ,1960
+ ,557
+ ,241
+ ,2103
+ ,999
+ ,612
+ ,1836
+ ,461
+ ,213
+ ,1864
+ ,561
+ ,264
+ ,1944
+ ,925
+ ,702
+ ,1935
+ ,471
+ ,297
+ ,1278
+ ,366
+ ,187
+ ,1744
+ ,660
+ ,292
+ ,2191
+ ,518
+ ,262
+ ,1893
+ ,598
+ ,274
+ ,2674
+ ,1526
+ ,1000
+ ,2617
+ ,307
+ ,203
+ ,2028
+ ,361
+ ,192
+ ,2412
+ ,745
+ ,465
+ ,2163
+ ,403
+ ,224
+ ,1920
+ ,404
+ ,316
+ ,2212
+ ,767
+ ,732
+ ,2319
+ ,565
+ ,347
+ ,1619
+ ,344
+ ,197
+ ,1746
+ ,571
+ ,344
+ ,2485
+ ,525
+ ,345
+ ,2079
+ ,557
+ ,361
+ ,2854
+ ,1604
+ ,1058
+ ,2651
+ ,374
+ ,236
+ ,2127
+ ,387
+ ,259
+ ,2154
+ ,644
+ ,404
+ ,2549
+ ,516
+ ,317
+ ,1912
+ ,443
+ ,287
+ ,2274
+ ,810
+ ,666
+ ,2197
+ ,533
+ ,434
+ ,1340
+ ,312
+ ,244
+ ,1952
+ ,560
+ ,404
+ ,2287
+ ,497
+ ,361
+ ,1667
+ ,475
+ ,342
+ ,2761
+ ,1445
+ ,1252
+ ,2092
+ ,332
+ ,254
+ ,1814
+ ,334
+ ,267
+ ,1919
+ ,750
+ ,552
+ ,1888
+ ,396
+ ,317
+ ,1514
+ ,413
+ ,352
+ ,1905
+ ,759
+ ,654
+ ,1870
+ ,493
+ ,455
+ ,1218
+ ,318
+ ,301
+ ,1830
+ ,612
+ ,439
+ ,2208
+ ,465
+ ,378
+ ,1759
+ ,455
+ ,404
+ ,2751
+ ,1485
+ ,1428
+ ,2455
+ ,327
+ ,326
+ ,1977
+ ,346
+ ,287
+ ,2512
+ ,705
+ ,662
+ ,2171
+ ,376
+ ,334
+ ,1772
+ ,390
+ ,316
+ ,2167
+ ,757
+ ,753
+ ,2237
+ ,469
+ ,443
+ ,1519
+ ,317
+ ,241
+ ,2023
+ ,580
+ ,442
+ ,2491
+ ,485
+ ,383
+ ,1881
+ ,456
+ ,445
+ ,3055
+ ,1566
+ ,1443
+ ,2653
+ ,328
+ ,272
+ ,2225
+ ,321
+ ,315
+ ,2462
+ ,682
+ ,687
+ ,2307
+ ,431
+ ,368
+ ,2186
+ ,430
+ ,451
+ ,2072
+ ,811
+ ,752
+ ,2151
+ ,455
+ ,462
+ ,1585
+ ,339
+ ,271
+ ,2092
+ ,592
+ ,553
+ ,2399
+ ,473
+ ,504
+ ,1882
+ ,458
+ ,497
+ ,2819
+ ,1891
+ ,1734
+ ,2267
+ ,278
+ ,292
+ ,1910
+ ,347
+ ,387
+ ,1975
+ ,652
+ ,727
+ ,1795
+ ,294
+ ,321
+ ,1549
+ ,393
+ ,429
+ ,1815
+ ,726
+ ,777
+ ,1742
+ ,472
+ ,549)
+ ,dim=c(3
+ ,91)
+ ,dimnames=list(c('OprichtingenVanVennootschappen'
+ ,'Kapitaalverhogingen'
+ ,'Kapitaalverminderingen')
+ ,1:91))
> y <- array(NA,dim=c(3,91),dimnames=list(c('OprichtingenVanVennootschappen','Kapitaalverhogingen','Kapitaalverminderingen'),1:91))
> 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 = '1'
> par3 <- 'No Linear Trend'
> par2 <- 'Do not include Seasonal 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, 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
OprichtingenVanVennootschappen Kapitaalverhogingen Kapitaalverminderingen
1 1925 358 155
2 1580 375 172
3 1961 761 467
4 1807 477 241
5 1526 547 294
6 1802 879 567
7 1822 450 280
8 1125 462 225
9 1569 1613 558
10 1829 854 342
11 1575 761 309
12 2339 1521 1437
13 2355 666 241
14 1960 557 241
15 2103 999 612
16 1836 461 213
17 1864 561 264
18 1944 925 702
19 1935 471 297
20 1278 366 187
21 1744 660 292
22 2191 518 262
23 1893 598 274
24 2674 1526 1000
25 2617 307 203
26 2028 361 192
27 2412 745 465
28 2163 403 224
29 1920 404 316
30 2212 767 732
31 2319 565 347
32 1619 344 197
33 1746 571 344
34 2485 525 345
35 2079 557 361
36 2854 1604 1058
37 2651 374 236
38 2127 387 259
39 2154 644 404
40 2549 516 317
41 1912 443 287
42 2274 810 666
43 2197 533 434
44 1340 312 244
45 1952 560 404
46 2287 497 361
47 1667 475 342
48 2761 1445 1252
49 2092 332 254
50 1814 334 267
51 1919 750 552
52 1888 396 317
53 1514 413 352
54 1905 759 654
55 1870 493 455
56 1218 318 301
57 1830 612 439
58 2208 465 378
59 1759 455 404
60 2751 1485 1428
61 2455 327 326
62 1977 346 287
63 2512 705 662
64 2171 376 334
65 1772 390 316
66 2167 757 753
67 2237 469 443
68 1519 317 241
69 2023 580 442
70 2491 485 383
71 1881 456 445
72 3055 1566 1443
73 2653 328 272
74 2225 321 315
75 2462 682 687
76 2307 431 368
77 2186 430 451
78 2072 811 752
79 2151 455 462
80 1585 339 271
81 2092 592 553
82 2399 473 504
83 1882 458 497
84 2819 1891 1734
85 2267 278 292
86 1910 347 387
87 1975 652 727
88 1795 294 321
89 1549 393 429
90 1815 726 777
91 1742 472 549
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Kapitaalverhogingen Kapitaalverminderingen
1746.72264 -0.09037 0.76321
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-751.70 -207.16 -24.08 216.19 757.96
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1746.72264 69.59055 25.100 < 2e-16 ***
Kapitaalverhogingen -0.09037 0.20709 -0.436 0.66364
Kapitaalverminderingen 0.76321 0.23616 3.232 0.00173 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 325.6 on 88 degrees of freedom
Multiple R-squared: 0.2869, Adjusted R-squared: 0.2707
F-statistic: 17.7 on 2 and 88 DF, p-value: 3.461e-07
> 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.2752955 0.5505911 0.7247045
[2,] 0.1413963 0.2827927 0.8586037
[3,] 0.4519836 0.9039672 0.5480164
[4,] 0.4096010 0.8192019 0.5903990
[5,] 0.3462452 0.6924905 0.6537548
[6,] 0.2718052 0.5436104 0.7281948
[7,] 0.1982356 0.3964713 0.8017644
[8,] 0.5765395 0.8469211 0.4234605
[9,] 0.5266968 0.9466064 0.4733032
[10,] 0.4765059 0.9530117 0.5234941
[11,] 0.3965445 0.7930891 0.6034555
[12,] 0.3254762 0.6509525 0.6745238
[13,] 0.2679852 0.5359704 0.7320148
[14,] 0.2121284 0.4242567 0.7878716
[15,] 0.3468442 0.6936885 0.6531558
[16,] 0.3048356 0.6096713 0.6951644
[17,] 0.3526071 0.7052142 0.6473929
[18,] 0.3055539 0.6111078 0.6944461
[19,] 0.4061807 0.8123615 0.5938193
[20,] 0.7557362 0.4885275 0.2442638
[21,] 0.7152322 0.5695355 0.2847678
[22,] 0.7398016 0.5203967 0.2601984
[23,] 0.7183758 0.5632484 0.2816242
[24,] 0.6605498 0.6789003 0.3394502
[25,] 0.5977744 0.8044512 0.4022256
[26,] 0.6013501 0.7972998 0.3986499
[27,] 0.5898373 0.8203253 0.4101627
[28,] 0.5742409 0.8515182 0.4257591
[29,] 0.6519487 0.6961026 0.3480513
[30,] 0.5962905 0.8074190 0.4037095
[31,] 0.6438677 0.7122645 0.3561323
[32,] 0.8263784 0.3472431 0.1736216
[33,] 0.7950744 0.4098512 0.2049256
[34,] 0.7524076 0.4951847 0.2475924
[35,] 0.8310023 0.3379954 0.1689977
[36,] 0.7904210 0.4191581 0.2095790
[37,] 0.7442342 0.5115316 0.2557658
[38,] 0.7017031 0.5965939 0.2982969
[39,] 0.8089849 0.3820301 0.1910151
[40,] 0.7676715 0.4646569 0.2323285
[41,] 0.7537464 0.4925071 0.2462536
[42,] 0.7538972 0.4922057 0.2461028
[43,] 0.7111710 0.5776581 0.2888290
[44,] 0.6700284 0.6599432 0.3299716
[45,] 0.6207822 0.7584356 0.3792178
[46,] 0.5921114 0.8157772 0.4078886
[47,] 0.5355275 0.9289449 0.4644725
[48,] 0.6086838 0.7826324 0.3913162
[49,] 0.6053388 0.7893223 0.3946612
[50,] 0.5647203 0.8705595 0.4352797
[51,] 0.7943505 0.4112990 0.2056495
[52,] 0.8090391 0.3819218 0.1909609
[53,] 0.7707918 0.4584163 0.2292082
[54,] 0.7701085 0.4597831 0.2298915
[55,] 0.7203465 0.5593070 0.2796535
[56,] 0.7870411 0.4259177 0.2129589
[57,] 0.7367238 0.5265525 0.2632762
[58,] 0.7294750 0.5410499 0.2705250
[59,] 0.6829317 0.6341367 0.3170683
[60,] 0.6721509 0.6556982 0.3278491
[61,] 0.6076457 0.7847087 0.3923543
[62,] 0.5553808 0.8892383 0.4446192
[63,] 0.6918278 0.6163444 0.3081722
[64,] 0.7209768 0.5580464 0.2790232
[65,] 0.6941021 0.6117958 0.3058979
[66,] 0.6564200 0.6871601 0.3435800
[67,] 0.6686631 0.6626738 0.3313369
[68,] 0.8150501 0.3698998 0.1849499
[69,] 0.7934049 0.4131901 0.2065951
[70,] 0.8176239 0.3647522 0.1823761
[71,] 0.8062451 0.3875099 0.1937549
[72,] 0.7840982 0.4318036 0.2159018
[73,] 0.7156965 0.5686070 0.2843035
[74,] 0.6645551 0.6708898 0.3354449
[75,] 0.7958810 0.4082381 0.2041190
[76,] 0.7234418 0.5531165 0.2765582
[77,] 0.8589701 0.2820597 0.1410299
[78,] 0.7646722 0.4706556 0.2353278
[79,] 0.6294249 0.7411502 0.3705751
[80,] 0.8243697 0.3512605 0.1756303
> postscript(file="/var/wessaorg/rcomp/tmp/10ea81353056152.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/24e8s1353056152.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/3jiia1353056152.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/4odve1353056152.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/5tvk51353056152.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 = 91
Frequency = 1
1 2 3 4 5 6
92.331070 -264.107242 -73.372330 -80.551268 -395.675673 -298.029937
7 8 9 10 11 12
-97.756227 -751.695433 -457.832765 -101.567408 -338.785559 -367.005707
13 14 15 16 17 18
484.527819 79.677975 -20.530406 -31.627309 -33.514333 -254.906122
19 20 21 22 23 24
4.166924 -578.368643 -165.937952 291.126363 -8.802883 301.967758
25 26 27 28 29 30
743.088472 167.363492 377.708237 281.736208 -31.388508 -24.080100
31 32 33 34 35 36
358.500914 -246.988759 -211.667271 522.412707 107.093086 444.750241
37 38 39 40 41 42
757.957118 217.578100 157.136969 606.969225 -13.731237 92.177307
43 44 45 46 47 48
167.210172 -564.751204 -52.453736 309.671154 -297.815948 189.319883
49 50 51 52 53 54
181.424032 -106.316933 -181.238980 -64.874640 -464.050685 -272.272846
55 56 57 58 59 60
-179.431805 -729.711833 -196.466987 214.804930 -254.942118 48.610000
61 62 63 64 65 66
489.021271 42.503305 323.741755 203.343523 -180.653626 -86.011111
67 68 69 70 71 72
194.557911 -383.009754 -8.648307 495.796204 -164.143256 348.481498
73 74 75 76 77 78
728.324837 266.874360 252.583163 318.364576 133.927996 -175.368164
79 80 81 82 83 84
92.791853 -337.917935 -23.279943 310.363721 -202.649310 -80.243058
85 86 87 88 89 90
322.542412 -100.727070 -267.656100 -170.144754 -489.624966 -459.129420
91
-381.070978
> postscript(file="/var/wessaorg/rcomp/tmp/6qgul1353056152.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 = 91
Frequency = 1
lag(myerror, k = 1) myerror
0 92.331070 NA
1 -264.107242 92.331070
2 -73.372330 -264.107242
3 -80.551268 -73.372330
4 -395.675673 -80.551268
5 -298.029937 -395.675673
6 -97.756227 -298.029937
7 -751.695433 -97.756227
8 -457.832765 -751.695433
9 -101.567408 -457.832765
10 -338.785559 -101.567408
11 -367.005707 -338.785559
12 484.527819 -367.005707
13 79.677975 484.527819
14 -20.530406 79.677975
15 -31.627309 -20.530406
16 -33.514333 -31.627309
17 -254.906122 -33.514333
18 4.166924 -254.906122
19 -578.368643 4.166924
20 -165.937952 -578.368643
21 291.126363 -165.937952
22 -8.802883 291.126363
23 301.967758 -8.802883
24 743.088472 301.967758
25 167.363492 743.088472
26 377.708237 167.363492
27 281.736208 377.708237
28 -31.388508 281.736208
29 -24.080100 -31.388508
30 358.500914 -24.080100
31 -246.988759 358.500914
32 -211.667271 -246.988759
33 522.412707 -211.667271
34 107.093086 522.412707
35 444.750241 107.093086
36 757.957118 444.750241
37 217.578100 757.957118
38 157.136969 217.578100
39 606.969225 157.136969
40 -13.731237 606.969225
41 92.177307 -13.731237
42 167.210172 92.177307
43 -564.751204 167.210172
44 -52.453736 -564.751204
45 309.671154 -52.453736
46 -297.815948 309.671154
47 189.319883 -297.815948
48 181.424032 189.319883
49 -106.316933 181.424032
50 -181.238980 -106.316933
51 -64.874640 -181.238980
52 -464.050685 -64.874640
53 -272.272846 -464.050685
54 -179.431805 -272.272846
55 -729.711833 -179.431805
56 -196.466987 -729.711833
57 214.804930 -196.466987
58 -254.942118 214.804930
59 48.610000 -254.942118
60 489.021271 48.610000
61 42.503305 489.021271
62 323.741755 42.503305
63 203.343523 323.741755
64 -180.653626 203.343523
65 -86.011111 -180.653626
66 194.557911 -86.011111
67 -383.009754 194.557911
68 -8.648307 -383.009754
69 495.796204 -8.648307
70 -164.143256 495.796204
71 348.481498 -164.143256
72 728.324837 348.481498
73 266.874360 728.324837
74 252.583163 266.874360
75 318.364576 252.583163
76 133.927996 318.364576
77 -175.368164 133.927996
78 92.791853 -175.368164
79 -337.917935 92.791853
80 -23.279943 -337.917935
81 310.363721 -23.279943
82 -202.649310 310.363721
83 -80.243058 -202.649310
84 322.542412 -80.243058
85 -100.727070 322.542412
86 -267.656100 -100.727070
87 -170.144754 -267.656100
88 -489.624966 -170.144754
89 -459.129420 -489.624966
90 -381.070978 -459.129420
91 NA -381.070978
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -264.107242 92.331070
[2,] -73.372330 -264.107242
[3,] -80.551268 -73.372330
[4,] -395.675673 -80.551268
[5,] -298.029937 -395.675673
[6,] -97.756227 -298.029937
[7,] -751.695433 -97.756227
[8,] -457.832765 -751.695433
[9,] -101.567408 -457.832765
[10,] -338.785559 -101.567408
[11,] -367.005707 -338.785559
[12,] 484.527819 -367.005707
[13,] 79.677975 484.527819
[14,] -20.530406 79.677975
[15,] -31.627309 -20.530406
[16,] -33.514333 -31.627309
[17,] -254.906122 -33.514333
[18,] 4.166924 -254.906122
[19,] -578.368643 4.166924
[20,] -165.937952 -578.368643
[21,] 291.126363 -165.937952
[22,] -8.802883 291.126363
[23,] 301.967758 -8.802883
[24,] 743.088472 301.967758
[25,] 167.363492 743.088472
[26,] 377.708237 167.363492
[27,] 281.736208 377.708237
[28,] -31.388508 281.736208
[29,] -24.080100 -31.388508
[30,] 358.500914 -24.080100
[31,] -246.988759 358.500914
[32,] -211.667271 -246.988759
[33,] 522.412707 -211.667271
[34,] 107.093086 522.412707
[35,] 444.750241 107.093086
[36,] 757.957118 444.750241
[37,] 217.578100 757.957118
[38,] 157.136969 217.578100
[39,] 606.969225 157.136969
[40,] -13.731237 606.969225
[41,] 92.177307 -13.731237
[42,] 167.210172 92.177307
[43,] -564.751204 167.210172
[44,] -52.453736 -564.751204
[45,] 309.671154 -52.453736
[46,] -297.815948 309.671154
[47,] 189.319883 -297.815948
[48,] 181.424032 189.319883
[49,] -106.316933 181.424032
[50,] -181.238980 -106.316933
[51,] -64.874640 -181.238980
[52,] -464.050685 -64.874640
[53,] -272.272846 -464.050685
[54,] -179.431805 -272.272846
[55,] -729.711833 -179.431805
[56,] -196.466987 -729.711833
[57,] 214.804930 -196.466987
[58,] -254.942118 214.804930
[59,] 48.610000 -254.942118
[60,] 489.021271 48.610000
[61,] 42.503305 489.021271
[62,] 323.741755 42.503305
[63,] 203.343523 323.741755
[64,] -180.653626 203.343523
[65,] -86.011111 -180.653626
[66,] 194.557911 -86.011111
[67,] -383.009754 194.557911
[68,] -8.648307 -383.009754
[69,] 495.796204 -8.648307
[70,] -164.143256 495.796204
[71,] 348.481498 -164.143256
[72,] 728.324837 348.481498
[73,] 266.874360 728.324837
[74,] 252.583163 266.874360
[75,] 318.364576 252.583163
[76,] 133.927996 318.364576
[77,] -175.368164 133.927996
[78,] 92.791853 -175.368164
[79,] -337.917935 92.791853
[80,] -23.279943 -337.917935
[81,] 310.363721 -23.279943
[82,] -202.649310 310.363721
[83,] -80.243058 -202.649310
[84,] 322.542412 -80.243058
[85,] -100.727070 322.542412
[86,] -267.656100 -100.727070
[87,] -170.144754 -267.656100
[88,] -489.624966 -170.144754
[89,] -459.129420 -489.624966
[90,] -381.070978 -459.129420
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -264.107242 92.331070
2 -73.372330 -264.107242
3 -80.551268 -73.372330
4 -395.675673 -80.551268
5 -298.029937 -395.675673
6 -97.756227 -298.029937
7 -751.695433 -97.756227
8 -457.832765 -751.695433
9 -101.567408 -457.832765
10 -338.785559 -101.567408
11 -367.005707 -338.785559
12 484.527819 -367.005707
13 79.677975 484.527819
14 -20.530406 79.677975
15 -31.627309 -20.530406
16 -33.514333 -31.627309
17 -254.906122 -33.514333
18 4.166924 -254.906122
19 -578.368643 4.166924
20 -165.937952 -578.368643
21 291.126363 -165.937952
22 -8.802883 291.126363
23 301.967758 -8.802883
24 743.088472 301.967758
25 167.363492 743.088472
26 377.708237 167.363492
27 281.736208 377.708237
28 -31.388508 281.736208
29 -24.080100 -31.388508
30 358.500914 -24.080100
31 -246.988759 358.500914
32 -211.667271 -246.988759
33 522.412707 -211.667271
34 107.093086 522.412707
35 444.750241 107.093086
36 757.957118 444.750241
37 217.578100 757.957118
38 157.136969 217.578100
39 606.969225 157.136969
40 -13.731237 606.969225
41 92.177307 -13.731237
42 167.210172 92.177307
43 -564.751204 167.210172
44 -52.453736 -564.751204
45 309.671154 -52.453736
46 -297.815948 309.671154
47 189.319883 -297.815948
48 181.424032 189.319883
49 -106.316933 181.424032
50 -181.238980 -106.316933
51 -64.874640 -181.238980
52 -464.050685 -64.874640
53 -272.272846 -464.050685
54 -179.431805 -272.272846
55 -729.711833 -179.431805
56 -196.466987 -729.711833
57 214.804930 -196.466987
58 -254.942118 214.804930
59 48.610000 -254.942118
60 489.021271 48.610000
61 42.503305 489.021271
62 323.741755 42.503305
63 203.343523 323.741755
64 -180.653626 203.343523
65 -86.011111 -180.653626
66 194.557911 -86.011111
67 -383.009754 194.557911
68 -8.648307 -383.009754
69 495.796204 -8.648307
70 -164.143256 495.796204
71 348.481498 -164.143256
72 728.324837 348.481498
73 266.874360 728.324837
74 252.583163 266.874360
75 318.364576 252.583163
76 133.927996 318.364576
77 -175.368164 133.927996
78 92.791853 -175.368164
79 -337.917935 92.791853
80 -23.279943 -337.917935
81 310.363721 -23.279943
82 -202.649310 310.363721
83 -80.243058 -202.649310
84 322.542412 -80.243058
85 -100.727070 322.542412
86 -267.656100 -100.727070
87 -170.144754 -267.656100
88 -489.624966 -170.144754
89 -459.129420 -489.624966
90 -381.070978 -459.129420
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/73jnz1353056152.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/83ng71353056152.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/wessaorg/rcomp/tmp/9l3g71353056152.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/wessaorg/rcomp/tmp/10tun31353056152.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/wessaorg/rcomp/tmp/113a6x1353056152.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/wessaorg/rcomp/tmp/12j2ks1353056152.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/wessaorg/rcomp/tmp/135md31353056152.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/wessaorg/rcomp/tmp/145rpp1353056152.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/wessaorg/rcomp/tmp/15j1e11353056152.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/wessaorg/rcomp/tmp/163u1v1353056152.tab")
+ }
>
> try(system("convert tmp/10ea81353056152.ps tmp/10ea81353056152.png",intern=TRUE))
character(0)
> try(system("convert tmp/24e8s1353056152.ps tmp/24e8s1353056152.png",intern=TRUE))
character(0)
> try(system("convert tmp/3jiia1353056152.ps tmp/3jiia1353056152.png",intern=TRUE))
character(0)
> try(system("convert tmp/4odve1353056152.ps tmp/4odve1353056152.png",intern=TRUE))
character(0)
> try(system("convert tmp/5tvk51353056152.ps tmp/5tvk51353056152.png",intern=TRUE))
character(0)
> try(system("convert tmp/6qgul1353056152.ps tmp/6qgul1353056152.png",intern=TRUE))
character(0)
> try(system("convert tmp/73jnz1353056152.ps tmp/73jnz1353056152.png",intern=TRUE))
character(0)
> try(system("convert tmp/83ng71353056152.ps tmp/83ng71353056152.png",intern=TRUE))
character(0)
> try(system("convert tmp/9l3g71353056152.ps tmp/9l3g71353056152.png",intern=TRUE))
character(0)
> try(system("convert tmp/10tun31353056152.ps tmp/10tun31353056152.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.583 1.155 7.784