t = [
2
3
4
5
6
7
10
15
20
30
];

spot = [
0.04859
0.04927
0.04952
0.05049
0.05152
0.05264
0.05323
0.05651
0.05777
0.0582
];

% QUESTION 2 & 3

% do a spline
ti = 0:.25:32;
spoti = spline (t, spot, ti);
g1 = plot (t, spot, 'o', ti, spoti);
% print -fg1 g1.ps


todaystr = '12-Jan-1998';
todaynum = datenum (todaystr);

format long
% first consider Six Canadian bond
Cmat = [
.0925, datenum('1-Dec-1999'), 107.65, .0493;
.0975, datenum('1-Dec-2001'), 116.42, .0503;
.075,  datenum('1-Dec-2003'), 111.89, .0513;
.0875, datenum('1-Dec-2005'), 122.53, .0522;
.07,   datenum('1-Dec-2006'), 112.26, .0525;
.095,  datenum('1-Jun-2010'), 136.46, .0542;
]

for i=1:6
  CanTheoPrice(i) = pvBond (todaystr, datestr (Cmat(i,2)), Cmat(i,1), 2, t, spot);
end
format short
CanTheoPrice

% now let's look at these Quebec bonds...
format long
Qmat = [
.1025, datenum('7-Apr-1998'),  101.05, 0.0538;
.06,   datenum('1-Apr-1999'),  101.34, 0.0484;
.1025, datenum('15-Oct-2001'), 117.05, 0.0518;
.0925, datenum('1-Apr-2002'),  114.95, 0.0524;
.075,  datenum('1-Dec-2003'),  110.33, 0.0542; 
.0937, datenum('16-Jan-2023'), 139.35, 0.0624;
]

for i=1:6
  QueTheoPrice(i) = pvBond (todaystr, datestr (Qmat(i,2)), Qmat(i,1), 2, t, spot);
end
format short
QueTheoPrice


% QUESTION 2
% first 3 col indicate which 3 data point to use, 
% last col is the Mean Squared Erro
M = [
1,2,3,  0;
1,2,4,  0;
1,2,5,  0;
1,2,6,  0;
1,2,7,  0;
1,2,8,  0;
1,2,9,  0;
1,2,10, 0;

1,3,4,  0;
1,3,5,  0;
1,3,6,  0;
1,3,7,  0;
1,3,8,  0;
1,3,9,  0;
1,3,10, 0;

1,4,5,  0;
1,4,6,  0;
1,4,7,  0;
1,4,8,  0;
1,4,9,  0;
1,4,10, 0;

1,5,6,  0;
1,5,7,  0;
1,5,8,  0;
1,5,9,  0;
1,5,10, 0;

1,6,7,  0;
1,6,8,  0;
1,6,9,  0;
1,6,10, 0;

1,7,8,  0;
1,7,9,  0;
1,7,10, 0;

1,8,9,  0;
1,8,10, 0;

1,9,10, 0;


2,3,4,  0;
2,3,5,  0;
2,3,6,  0;
2,3,7,  0;
2,3,8,  0;
2,3,9,  0;
2,3,10, 0;

2,4,5,  0;
2,4,6,  0;
2,4,7,  0;
2,4,8,  0;
2,4,9,  0;
2,4,10, 0;

2,5,6,  0;
2,5,7,  0;
2,5,8,  0;
2,5,9,  0;
2,5,10, 0;

2,6,7,  0;
2,6,8,  0;
2,6,9,  0;
2,6,10, 0;

2,7,8,  0;
2,7,9,  0;
2,7,10, 0;

2,8,9,  0;
2,8,10, 0;

2,9,10, 0;

3,4,5,  0;
3,4,6,  0;
3,4,7,  0;
3,4,8,  0;
3,4,9,  0;
3,4,10, 0;

3,5,6,  0;
3,5,7,  0;
3,5,8,  0;
3,5,9,  0;
3,5,10, 0;

3,6,7,  0;
3,6,8,  0;
3,6,9,  0;
3,6,10, 0;

3,7,8,  0;
3,7,9,  0;
3,7,10, 0;

3,8,9,  0;
3,8,10, 0;

3,9,10, 0;


4,5,6,  0;
4,5,7,  0;
4,5,8,  0;
4,5,9,  0;
4,5,10, 0;

4,6,7,  0;
4,6,8,  0;
4,6,9,  0;
4,6,10, 0;

4,7,8,  0;
4,7,9,  0;
4,7,10, 0;

4,8,9,  0;
4,8,10, 0;

4,9,10, 0;

5,6,7,  0;
5,6,8,  0;
5,6,9,  0;
5,6,10, 0;

5,7,8,  0;
5,7,9,  0;
5,7,10, 0;

5,8,8,  0;
5,8,9,  0;

5,9,10, 0;


6,7,8,  0;
6,7,9,  0;
6,7,10, 0;

6,8,9,  0;
6,8,10, 0;

6,9,10, 0;


7,8,9,  0;
7,8,10, 0;

7,9,10, 0;

8,9,10, 0;

];

pos = -10;
min_mse = 100000;
for i = 1:length(M(:,1))
  bb = mySolve (M,i, t, spot);
  M(i,4) = mse(M, t, spot, bb);

  if min_mse > M(i,4)
     min_mse = M(i,4);
     pos = i;
     b = bb;
  end
end

format long
b
M(pos,:)


% consider first bond