% LMCS010201_01.m
% Get Standard & Poors stock prices 2/17/94 to 2/9/96 from 
% http://finance.yahoo.com/q?s=^GSPC&d=c&t=6m&l=on&z=m&q=l
% data used to generate an X, Y plot for Homework
% See lmcs010201_2010.m for 2010 data
% See lmcs010201_01.m for 2001 data
% These data copied from csv file from most recent to oldest
% must be flipped for runs test.
% Note that the book just uses runs of downturn in analysis.
% None of Gallagher's programs are written to handle downturns, so
% several m.files will have to be rewritten to match the book.
% Revised 9/15/10 for EEOS601
% Written by Eugene.Gallagher@umb.edu for EEOS601
SandP500=[656.37
656.07
649.93
646.33
641.43
635.84
638.46
636.02
630.15
624.22
621.62
617.03
619.96
612.79
613.4
611.83
608.24
606.37
608.44
599.82
601.81
602.69
598.48
609.45
618.46
616.71
617.7
621.32
620.73
615.93
614.12
614.53
614.3
611.95
610.49
605.94
611.93
606.81
616.34
616.92
621.69
618.78
619.52
617.48
616.17
620.18
617.68
613.68
606.98
605.37
607.64
606.45
601.32
599.97
598.4
600.24
596.85
600.07
597.34
593.96
589.29
592.3
592.72
593.26
591.71
586.32
588.46
590.57
589.72
584.22
581.5
583.25
579.7
576.72
582.47
586.54
585.06
587.46
590.65
587.44
586.78
583.03
584.5
583.1
579.46
577.52
578.37
582.49
582.63
581.47
582.34
581.72
584.41
585.87
581.04
581.41
581.81
581.73
583
586.77
584.2
582.77
583.35
583.61
578.77
576.51
573.91
572.68
570.29
570.17
569.17
563.84
561.88
560.92
560
559.05
560.1
557.46
557.14
559.52
558.11
559.21
559.04
559.97
558.57
559.74
555.11
557.45
559.71
560.39
560.03
558.94
558.75
558.8
559.64
562.06
562.93
565.22
561.61
561.1
556.63
553.62
553.54
550.98
558.46
562.72
559.89
561
560.89
554.78
557.19
556.37
553.99
547.26
547.09
544.75
543.87
544.73
542.43
544.13
549.71
551.07
543.98
544.98
545.22
539.83
537.12
536.47
536.05
530.88
527.94
532.35
533.13
535.55
535.6
532.51
533.49
533.4
523.58
523.65
528.59
528.61
528.59
523.65
519.19
519.58
527.07
528.19
527.74
525.55
524.37
524.36
523.56
523.96
520.12
520.54
520.48
514.86
514.26
514.71
513.55
512.66
512.1
512.89
508.49
505.29
504.92
505.37
506.13
509.23
507.17
505.53
507.01
506.42
506.08
505.57
505.24
501.85
500.71
502.22
503.12
503.9
503.2
500.97
495.95
495.67
495.07
496.14
495.52
495.41
491.88
492.89
490.05
489.57
483.16
483.14
482.12
485.63
485.42
485.13
485.65
487.39
483.81
488.11
486.91
485.07
482.72
481.97
485.22
484.54
482.55
481.65
481.46
480.19
481.19
480.81
481.14
478.65
472.79
470.4
470.42
468.51
470.39
468.32
467.44
465.86
465.82
464.78
466.95
469.71
470.05
469.38
465.97
461.64
461.66
461.68
460.83
460.68
460.34
460.71
459.11
459.27
461.17
460.86
462.47
459.83
459.68
459.61
457.1
457.91
458.8
455.34
454.97
450.15
449.47
446.96
445.45
451.23
453.11
453.32
453.3
448.92
453.69
455.17
454.16
452.29
449.93
450.09
458.3
461.47
463.57
465.62
465.03
466.04
462.35
464.37
465.4
465.65
463.07
462.28
467.91
466.51
468.42
472.35
473.77
465.85
462.62
461.53
460.83
464.89
466.85
470.28
467.66
468.96
469.1
467.79
465.47
465.79
459.04
455.1
452.36
453.52
454.59
461.74
462.71
462.24
464.84
462.05
460.82
459.67
461.27
461.46
463.36
470.85
471.19
474.81
468.8
467.51
466.21
468.18
473.14
470.99
471.86
470.99
473.17
475.49
476.07
474.59
473.8
468.08
469.03
464.51
462.32
463.68
463.17
465.17
465.01
461.23
461.94
458.88
460.3
457.92
457.89
457.09
458.4
461.45
460.56
461.01
458.26
454.24
452.57
453.36
454.25
453.11
452.61
451.6
453.86
455.22
454.16
453.41
448.73
447.95
448.06
449.55
448.38
446.13
446.37
446.2
444.27
447.63
446.07
447.31
442.8
449.63
453.09
451.34
455.48
458.45
461.93
460.61
462.37
459.1
458.67
457.86
457.06
458.21
458.88
460.13
457.65
457.63
456.5
457.33
457.06
456.34
454.81
453.2
454.92
456.48
453.69
449.37
444.49
444.14
443.75
441.49
446.01
442.32
447.82
451.38
451.72
453.03
453.02
450.91
449.1
451.87
452.71
447.63
448.73
441.96
442.54
442.46
446.18
446.38
446.26
447.57
449.87
447.1
450.88
448.05
448.29
438.92
445.77
445.55
452.48
460
460.58
464.35
468.54
468.8
468.54
471.06
470.9
469.42
467.01
467.39
466.44
463.9
467.06
465.88
466.91
464.74
463.01
464.81
464.44
467.14
466.07
464.26
470.69
471.46
467.69
470.34]
% cday is calendar day with day 1 being 2/17/94, dates run from 1 to 723
cday =[723
722
721
720
719
716
715
714
713
712
709
708
707
706
705
702
701
700
699
698
695
694
693
692
691
688
687
686
685
681
680
679
678
674
673
672
671
670
667
666
665
664
663
660
659
658
657
656
653
652
651
650
649
646
644
643
642
639
638
637
636
635
632
631
630
629
628
625
624
623
622
621
618
617
616
615
614
611
610
609
608
607
604
603
602
601
600
597
596
595
594
593
590
589
588
587
586
583
582
581
580
579
576
575
574
573
572
569
568
567
566
562
561
560
559
558
555
554
553
552
551
548
547
546
545
544
541
540
539
538
537
534
533
532
531
530
527
526
525
524
523
520
519
518
517
516
513
512
511
510
509
506
505
504
502
499
498
497
496
495
492
491
490
489
488
485
484
483
482
481
478
477
476
475
474
471
470
469
468
464
463
462
461
460
457
456
455
454
453
450
449
448
447
446
443
442
441
440
439
436
435
434
433
432
429
428
427
426
425
421
420
419
418
415
414
413
412
411
408
407
406
405
404
401
400
399
398
397
394
393
392
391
390
387
386
385
384
383
380
379
378
377
376
373
372
371
370
366
365
364
363
362
359
358
357
356
355
352
351
350
349
348
345
344
343
342
341
338
337
336
335
334
331
330
329
328
327
324
323
322
321
317
316
315
314
310
309
308
307
306
303
302
301
300
299
296
295
294
293
292
289
288
287
286
285
282
280
279
278
275
274
273
272
271
268
267
266
265
264
261
260
259
258
257
254
253
252
251
250
247
246
245
244
243
240
239
238
237
236
233
232
231
230
229
226
225
224
223
222
219
218
217
216
215
212
211
210
209
208
205
204
203
202
198
197
196
195
194
191
190
189
188
187
184
183
182
181
180
177
176
175
174
173
170
169
168
167
166
163
162
161
160
159
156
155
154
153
152
149
148
147
146
145
142
141
140
139
135
134
133
132
131
128
127
126
125
124
121
120
119
118
117
114
113
112
111
110
107
106
105
104
100
99
98
97
96
93
92
91
90
89
86
85
84
83
82
79
78
77
76
75
72
71
69
68
65
64
63
62
61
58
57
56
55
54
51
50
49
48
47
43
42
41
40
37
36
35
34
33
30
29
28
27
26
23
22
21
20
19
16
15
14
13
12
9
8
7
6
2
1]
whos
SandP=flipud([cday SandP500])
[r,c]=size(SandP);
tradingday=[1:r]'
SandP=flipud([tradingday cday SandP500])
% Create character matrix of calendar dates using built-in functions.
callabel=datestr(SandP(:,2)+datenum('16.02.1994','dd.mm.yyyy'),2)

% Use Gallagher's runs2 to analyze the observed number of runs
[obsruns,out]=runs2(SandP500)
% obsruns contains the observed number of runs

% Using the formula on page 12 (L&M 4th edition), calculate the expected
% number of runs if the probability of the S & P going ups
% is 0.5:
[r,c]=size(obsruns);
k=[1:r]';   % Solve for runs of length 1, 2, 3, ... max run length.
p=.5;  % probability of the run going up (or down)
T=out(1); % Total number of runs
% Using the formula on page 13, calculate the expected number of runs if
% the probability of the S & P going ups
% is 0.5:

Expected=T*p.^(k-1)*(1-p)  
% This is the formula immediately above Table 
% 1.2.2 on page 13 in Larsen & Marx.

% Comparing the Observed and Expected, they seem quite close, with maybe too many 2-d runs:
disp('No. of Runs Observed  Expected')
disp([obsruns Expected])

% Goodness of fit tests. To perform the chi-square test, we need to pool 
% the categories so that each has an expected value of 5 or more runs.
E=Expected(1:5);  % Expected Number of runs of lengths 1, 2, 3, 4
E(6)=sum(Expected(6)+Expected(7));    % The expected number of runs of length 6 7 7
                  % since we are modeling a problem with 242 total runs.
E(7)=242-sum(E);
O=[obsruns(1:5,2);obsruns(6,2)+obsruns(7,2)];   % Observed runs
O(7)=242-sum(O);    % The number of observed runs of 5 or greater can be summed or obtained by difference.
disp([O E])

% Chi-square statistic, sum of observed minus expected squared, divided by expected.

Chi2=sum ((O-E).^2 ./ E)

df=6    % Degrees of freedom for this problem is no. categories minus 1
% without the statistics toolbox use this: p = 1-gammainc(Chi2/2, df/2)
% With the statistics toolbox use this cumulative distribution function
p = 1-chi2cdf(Chi2,6)

% The p value, is 0.1145, which indicates there is only weak evidence 
% against the random model
% The probability of observing a result as extreme or more extreme than the
% one observed is about 11%
% By convention (dating back to Fisher), only p values less than about 0.05
% are considered moderate to strong evidence that the null hypothesis is
% false.
% Note that the book's Table 1.2.2 presented what is the 'best' estimate
% for p. Can we get a better fit to the data by choosing a slightly 
% different value of p?  The best value of p will have the lowest chi square
% statistic, so let's find how chi-square varies as a function of p.
pest=0.52:0.00001:0.522; % We'll check all of the p values from 0.2 to 0.8
lp=length(pest);    % i gives us the number of values that we'll be testing
Chi2est=zeros(lp,1); % Create a matrix of all zeros to save the results of
% the simulation.
O=O(1:6);O(6)=242-sum(O(1:5));
for i=1:lp
    k=[1:12]';   % Solve for runs of length 1, 2, 3, ... 12 simultaneously.
    p=pest(i);  % insert the probability from the pest vector
    T=242;       % This doesn't need to be inside the for loop since its
                 % value is unchanged, but it doesn't hurt.
    E=T*p.^(k-1)*(1-p);  % Calculate the Number of runs of length 1 to 12 
    % using equation on p. 13 in Larsen & Marx
    E=E(1:6);
    E(6)=T-sum(E(1:5)); % Get the probability of 6 or more runs by difference
    Chi2est(i)=sum ((O-E).^2 ./ E);
end
plot(pest,Chi2est);grid;xlabel('Estimated p');ylabel('Chi-square');figure(gcf);pause
% From this plot, it appears that a p value of about 0.5214 would provide a better fit than 0.5
% We can find the nearly exact value of p that produces a minimum of chi square using ginput
% this can also be solved numerically using the Matlab function fzero.
 [x,y]=ginput(1)
 % When uncommented ginput puts up cross hairs to find the minimum, which is about 0.46
% Analysis of the model with p=0.5214
    k=[1:12]';   % Solve for runs of length 1, 2, 3, 4 simultaneously.
    p=0.5212;  % insert the probability from the pest vector
    T=242;       % This doesn't need to be inside the for loop since its value is unchanged, but it doesn't hurt.
    E=T*p.^(k-1)*(1-p);  % Calculate the Number of runs of length 1 to 4 using equation on p. 13 in Larsen & Marx
    E=E(1:6);E(6)=242-sum(E(1:5));    % Get the probability of 5 or more runs by difference
    disp('   Observed  Expected (given p=0.5214)')
    disp([O E])
    Chi2est=sum ((O-E).^2 ./ E)
    df=5
    p=1-gammainc(Chi2est/2,df/2)
    % this produces a p value for the chi-square test of 0.1554, which is a little bit better fit to the model
    % The interpretation of this is that if the S & P went up one day, these data indicate that there is only a 46%
    % chance that it will go up the second day (Compare Table 1.2.1 to see that for p values less than 0.5, the
    % probability of short runs is enhanced and long runs reduced)
