# linear Gaussian state space model - local level model
# clear the environment
rm(list=ls())
# annual flow of the Nile from 1871-1970
data(Nile)
help(Nile)
plot(Nile)
y <- datasets::Nile
# uncomment the following two line to see the analysis with missing values
# y[21:50] <- NA
# y[71:80] <- NA
# sample size
T <-length(y)
library(KFAS)
## Warning: package 'KFAS' was built under R version 3.2.4
# define the state-space local level model
y.LLM <- SSModel(y ~ SSMtrend(1, Q=list(NA)), H=NA )
# check state variable
y.LLM$a
## [,1]
## level 0
# check system matrices
y.LLM$Z
## , , 1
##
## level
## [1,] 1
y.LLM$T
## , , 1
##
## level
## level 1
y.LLM$R
## , , 1
##
## [,1]
## level 1
y.LLM$H
## , , 1
##
## [,1]
## [1,] NA
y.LLM$Q
## , , 1
##
## [,1]
## [1,] NA
# maximum likelihood estimation of paramaters of Q and H (i.e. variance of the two inovations)
pars <- log( rep(var(y, na.rm=TRUE), 2) )
y.fit <- fitSSM( model=y.LLM, inits=pars, method='BFGS' )
y.fit
## $optim.out
## $optim.out$par
## [1] 7.292448 9.622361
##
## $optim.out$value
## [1] 632.5456
##
## $optim.out$counts
## function gradient
## 35 10
##
## $optim.out$convergence
## [1] 0
##
## $optim.out$message
## NULL
##
##
## $model
## Call:
## SSModel(formula = y ~ SSMtrend(1, Q = list(NA)), H = NA)
##
## State space model object of class SSModel
##
## Dimensions:
## [1] Number of time points: 100
## [1] Number of time series: 1
## [1] Number of disturbances: 1
## [1] Number of states: 1
## Names of the states:
## [1] level
## Distributions of the time series:
## [1] gaussian
##
## Object is a valid object of class SSModel.
# run Kalman Filter and Smoother with estimated parameters
y.out <- KFS(y.fit$model)
# construct 90% confidence intervals for filtered state - shorter way, using predict function with missing n.ahead option
y.KF.CI <- predict(y.fit$model, interval="confidence", level=0.9, filtered=TRUE)
# construct 90% confidence intervals for smoothed state
y.KS.CI <- predict(y.fit$model, interval="confidence", level=0.9)
# construct 90% confidence intervals for filtered state - longer, do it yourself way
y.KF.CI_2 <- as.vector(y.out$a) + sqrt(cbind(y.out$P)) %*% as.vector( qnorm(0.95)*c(-1,1) )
# construct 90% confidence intervals for smoothed state
y.KS.CI_2 <- as.vector(y.out$alphahat) + sqrt(cbind(y.out$V)) %*% as.vector( qnorm(0.95)*c(-1,1) )
# replace the filtered state for first period by NA
y.KF.CI[1,] <- NA
# check that confidence intervals are actually same for both methods
cbind( y.KF.CI, y.KF.CI_2[-(T+1),] )
## Time Series:
## Start = 1871
## End = 1970
## Frequency = 1
## y.KF.CI.fit y.KF.CI.lwr y.KF.CI.upr y.KF.CI_2[-(T + 1), ].Series 1
## 1871 NA NA NA 0.0000
## 1872 1120.0000 908.2810 1331.7190 908.2810
## 1873 1140.9279 981.7188 1300.1370 981.7188
## 1874 1072.7980 932.7385 1212.8575 932.7385
## 1875 1117.3092 986.0556 1248.5628 986.0556
## 1876 1129.9725 1003.1175 1256.8276 1003.1175
## 1877 1138.4585 1013.8818 1263.0351 1013.8818
## 1878 1048.8549 925.4794 1172.2304 925.4794
## 1879 1098.0296 975.2930 1220.7662 975.2930
## 1880 1171.3031 1048.9078 1293.6983 1048.9078
## 1881 1162.9038 1040.6914 1285.1161 1040.6914
## 1882 1117.9500 995.8358 1240.0643 995.8358
## 1883 1069.0254 946.9639 1191.0870 946.9639
## 1884 1079.9760 957.9427 1202.0093 957.9427
## 1885 1057.0066 934.9884 1179.0247 934.9884
## 1886 1047.1217 925.1117 1169.1316 925.1117
## 1887 1023.8527 901.8471 1145.8583 901.8471
## 1888 1065.5552 943.5520 1187.5585 943.5520
## 1889 994.3679 872.3659 1116.3698 872.3659
## 1890 984.6555 862.6542 1106.6567 862.6542
## 1891 1026.1415 904.1406 1148.1424 904.1406
## 1892 1045.8659 923.8652 1167.8666 923.8652
## 1893 1089.6989 967.6983 1211.6995 967.6983
## 1894 1105.8027 983.8021 1227.8032 983.8021
## 1895 1144.3114 1022.3109 1266.3119 1022.3109
## 1896 1175.2067 1053.2062 1297.2072 1053.2062
## 1897 1187.1690 1065.1685 1309.1695 1065.1685
## 1898 1145.1961 1023.1957 1267.1966 1023.1957
## 1899 1133.1263 1011.1258 1255.1268 1011.1258
## 1900 1037.2196 915.2191 1159.2201 915.2191
## 1901 984.5511 862.5506 1106.5515 862.5506
## 1902 955.0278 833.0273 1077.0283 833.0273
## 1903 885.3189 763.3184 1007.3194 763.3184
## 1904 899.9218 777.9213 1021.9223 777.9213
## 1905 882.0500 760.0495 1004.0504 760.0495
## 1906 833.6996 711.6991 955.7001 711.6991
## 1907 855.6784 733.6779 977.6788 733.6779
## 1908 811.9672 689.9667 933.9676 689.9667
## 1909 867.5235 745.5230 989.5239 745.5230
## 1910 916.2548 794.2543 1038.2553 794.2543
## 1911 930.3407 808.3402 1052.3412 808.3402
## 1912 903.8112 781.8107 1025.8117 781.8107
## 1913 856.3258 734.3253 978.3263 734.3253
## 1914 749.4166 627.4161 871.4171 627.4161
## 1915 769.3345 647.3340 891.3350 647.3340
## 1916 751.3524 629.3520 873.3529 629.3520
## 1917 849.8018 727.8013 971.8023 727.8013
## 1918 916.6186 794.6181 1038.6191 794.6181
## 1919 894.0207 772.0202 1016.0212 772.0202
## 1920 859.2980 737.2975 981.2985 737.2975
## 1921 849.0703 727.0698 971.0708 727.0698
## 1922 827.4200 705.4195 949.4205 705.4195
## 1923 832.1149 710.1144 954.1153 710.1144
## 1924 840.6300 718.6295 962.6304 718.6295
## 1925 846.3369 724.3365 968.3374 724.3365
## 1926 806.7228 684.7223 928.7232 684.7223
## 1927 816.9449 694.9444 938.9454 694.9444
## 1928 797.4646 675.4641 919.4651 675.4641
## 1929 797.0734 675.0730 919.0739 675.0730
## 1930 861.9483 739.9478 983.9488 739.9478
## 1931 834.4554 712.4549 956.4559 712.4549
## 1932 820.1798 698.1794 942.1803 698.1794
## 1933 832.1493 710.1488 954.1498 710.1488
## 1934 835.5812 713.5807 957.5816 713.5807
## 1935 864.5350 742.5345 986.5355 742.5345
## 1936 896.4388 774.4383 1018.4393 774.4383
## 1937 896.5887 774.5882 1018.5891 774.5882
## 1938 876.6693 754.6689 998.6698 754.6689
## 1939 912.2760 790.2755 1034.2765 790.2755
## 1940 874.5475 752.5470 996.5480 752.5470
## 1941 821.5243 699.5238 943.5248 699.5238
## 1942 775.4507 653.4503 897.4512 653.4503
## 1943 794.2913 672.2908 916.2918 672.2908
## 1944 799.0205 677.0200 921.0210 677.0200
## 1945 783.7929 661.7924 905.7934 661.7924
## 1946 788.3881 666.3876 910.3886 666.3876
## 1947 855.5825 733.5820 977.5829 733.5820
## 1948 856.7622 734.7617 978.7627 734.7617
## 1949 861.3656 739.3652 983.3661 739.3652
## 1950 857.7963 735.7958 979.7968 735.7958
## 1951 866.3965 744.3960 988.3969 744.3960
## 1952 833.7098 711.7093 955.7103 711.7093
## 1953 811.0876 689.0871 933.0881 689.0871
## 1954 818.2747 696.2742 940.2752 696.2742
## 1955 880.1582 758.1577 1002.1587 758.1577
## 1956 890.2641 768.2636 1012.2645 768.2636
## 1957 915.8309 793.8304 1037.8314 793.8304
## 1958 884.0964 762.0960 1006.0969 762.0960
## 1959 894.4858 772.4854 1016.4863 772.4854
## 1960 915.9876 793.9871 1037.9881 793.9871
## 1961 889.0183 767.0178 1011.0188 767.0178
## 1962 923.9977 801.9972 1045.9982 801.9972
## 1963 919.1913 797.1908 1041.1918 797.1908
## 1964 914.3332 792.3327 1036.3337 792.3327
## 1965 982.6104 860.6099 1104.6109 860.6099
## 1966 963.7535 841.7530 1085.7540 841.7530
## 1967 905.6013 783.6008 1027.6017 783.6008
## 1968 909.1795 787.1790 1031.1799 787.1790
## 1969 858.1239 736.1235 980.1244 736.1235
## 1970 819.6349 697.6344 941.6353 697.6344
## y.KF.CI_2[-(T + 1), ].Series 2
## 1871 0.0000
## 1872 1331.7190
## 1873 1300.1370
## 1874 1212.8575
## 1875 1248.5628
## 1876 1256.8276
## 1877 1263.0351
## 1878 1172.2304
## 1879 1220.7662
## 1880 1293.6983
## 1881 1285.1161
## 1882 1240.0643
## 1883 1191.0870
## 1884 1202.0093
## 1885 1179.0247
## 1886 1169.1316
## 1887 1145.8583
## 1888 1187.5585
## 1889 1116.3698
## 1890 1106.6567
## 1891 1148.1424
## 1892 1167.8666
## 1893 1211.6995
## 1894 1227.8032
## 1895 1266.3119
## 1896 1297.2072
## 1897 1309.1695
## 1898 1267.1966
## 1899 1255.1268
## 1900 1159.2201
## 1901 1106.5515
## 1902 1077.0283
## 1903 1007.3194
## 1904 1021.9223
## 1905 1004.0504
## 1906 955.7001
## 1907 977.6788
## 1908 933.9676
## 1909 989.5239
## 1910 1038.2553
## 1911 1052.3412
## 1912 1025.8117
## 1913 978.3263
## 1914 871.4171
## 1915 891.3350
## 1916 873.3529
## 1917 971.8023
## 1918 1038.6191
## 1919 1016.0212
## 1920 981.2985
## 1921 971.0708
## 1922 949.4205
## 1923 954.1153
## 1924 962.6304
## 1925 968.3374
## 1926 928.7232
## 1927 938.9454
## 1928 919.4651
## 1929 919.0739
## 1930 983.9488
## 1931 956.4559
## 1932 942.1803
## 1933 954.1498
## 1934 957.5816
## 1935 986.5355
## 1936 1018.4393
## 1937 1018.5891
## 1938 998.6698
## 1939 1034.2765
## 1940 996.5480
## 1941 943.5248
## 1942 897.4512
## 1943 916.2918
## 1944 921.0210
## 1945 905.7934
## 1946 910.3886
## 1947 977.5829
## 1948 978.7627
## 1949 983.3661
## 1950 979.7968
## 1951 988.3969
## 1952 955.7103
## 1953 933.0881
## 1954 940.2752
## 1955 1002.1587
## 1956 1012.2645
## 1957 1037.8314
## 1958 1006.0969
## 1959 1016.4863
## 1960 1037.9881
## 1961 1011.0188
## 1962 1045.9982
## 1963 1041.1918
## 1964 1036.3337
## 1965 1104.6109
## 1966 1085.7540
## 1967 1027.6017
## 1968 1031.1799
## 1969 980.1244
## 1970 941.6353
cbind( y.KS.CI, y.KS.CI )
## Time Series:
## Start = 1871
## End = 1970
## Frequency = 1
## y.KS.CI.fit y.KS.CI.lwr y.KS.CI.upr y.KS.CI.fit y.KS.CI.lwr
## 1871 1111.6686 1007.2213 1216.1159 1111.6686 1007.2213
## 1872 1110.8579 1017.1889 1204.5270 1110.8579 1017.1889
## 1873 1105.2655 1017.9341 1192.5969 1105.2655 1017.9341
## 1874 1113.5161 1029.7872 1197.2451 1113.5161 1029.7872
## 1875 1112.3785 1030.6504 1194.1066 1112.3785 1030.6504
## 1876 1106.6070 1025.9743 1187.2398 1106.6070 1025.9743
## 1877 1095.6402 1015.6021 1175.6783 1095.6402 1015.6021
## 1878 1112.1755 1032.4587 1191.8924 1112.1755 1032.4587
## 1879 1117.2460 1037.7022 1196.7897 1117.2460 1037.7022
## 1880 1097.7224 1018.2718 1177.1730 1097.7224 1018.2718
## 1881 1074.0850 994.6845 1153.4855 1074.0850 994.6845
## 1882 1058.1430 978.7694 1137.5165 1058.1430 978.7694
## 1883 1054.1832 974.8241 1133.5423 1054.1832 974.8241
## 1884 1044.7923 965.4410 1124.1436 1044.7923 965.4410
## 1885 1040.3437 960.9965 1119.6908 1040.3437 960.9965
## 1886 1037.8746 958.5297 1117.2195 1037.8746 958.5297
## 1887 1042.9830 963.6393 1122.3267 1042.9830 963.6393
## 1888 1034.7591 955.4161 1114.1022 1034.7591 955.4161
## 1889 1049.4756 970.1329 1128.8183 1049.4756 970.1329
## 1890 1073.0930 993.7505 1152.4355 1073.0930 993.7505
## 1891 1090.2001 1010.8577 1169.5425 1090.2001 1010.8577
## 1892 1106.3536 1027.0112 1185.6959 1106.3536 1027.0112
## 1893 1112.4219 1033.0795 1191.7642 1112.4219 1033.0795
## 1894 1114.8336 1035.4913 1194.1759 1114.8336 1035.4913
## 1895 1104.0931 1024.7508 1183.4354 1104.0931 1024.7508
## 1896 1078.1822 998.8399 1157.5245 1078.1822 998.8399
## 1897 1038.4718 959.1295 1117.8141 1038.4718 959.1295
## 1898 999.5858 920.2435 1078.9281 999.5858 920.2435
## 1899 950.9290 871.5867 1030.2713 950.9290 871.5867
## 1900 919.4882 840.1459 998.8305 919.4882 840.1459
## 1901 895.7819 816.4396 975.1242 895.7819 816.4396
## 1902 874.1950 794.8527 953.5373 874.1950 794.8527
## 1903 870.1419 790.7996 949.4842 870.1419 790.7996
## 1904 859.2914 779.9491 938.6337 859.2914 779.9491
## 1905 850.9990 771.6567 930.3413 850.9990 771.6567
## 1906 857.3023 777.9600 936.6446 857.3023 777.9600
## 1907 857.8940 778.5517 937.2363 857.8940 778.5517
## 1908 874.6278 795.2855 953.9701 874.6278 795.2855
## 1909 877.2164 797.8741 956.5587 877.2164 797.8741
## 1910 862.9923 783.6500 942.3346 862.9923 783.6500
## 1911 838.4533 759.1110 917.7956 838.4533 759.1110
## 1912 814.6395 735.2972 893.9818 814.6395 735.2972
## 1913 799.4507 720.1084 878.7930 799.4507 720.1084
## 1914 817.6811 738.3388 897.0234 817.6811 738.3388
## 1915 835.2967 755.9544 914.6390 835.2967 755.9544
## 1916 865.8826 786.5403 945.2249 865.8826 786.5403
## 1917 871.7418 792.3995 951.0841 871.7418 792.3995
## 1918 855.3905 776.0482 934.7328 855.3905 776.0482
## 1919 841.3152 761.9729 920.6575 841.3152 761.9729
## 1920 834.7630 755.4207 914.1053 834.7630 755.4207
## 1921 829.5500 750.2077 908.8923 829.5500 750.2077
## 1922 830.3261 750.9838 909.6684 830.3261 750.9838
## 1923 829.6744 750.3321 909.0166 829.6744 750.3321
## 1924 825.6826 746.3403 905.0249 825.6826 746.3403
## 1925 818.1569 738.8146 897.4992 818.1569 738.8146
## 1926 822.3231 742.9808 901.6654 822.3231 742.9808
## 1927 824.2827 744.9404 903.6250 824.2827 744.9404
## 1928 834.0541 754.7119 913.3964 834.0541 754.7119
## 1929 847.5284 768.1861 926.8707 847.5284 768.1861
## 1930 842.2744 762.9321 921.6167 842.2744 762.9321
## 1931 845.1233 765.7810 924.4656 845.1233 765.7810
## 1932 854.2117 774.8694 933.5540 854.2117 774.8694
## 1933 862.2504 782.9081 941.5927 862.2504 782.9081
## 1934 871.9676 792.6253 951.3099 871.9676 792.6253
## 1935 874.6757 795.3334 954.0180 874.6757 795.3334
## 1936 866.7461 787.4038 946.0884 866.7461 787.4038
## 1937 855.8727 776.5304 935.2150 855.8727 776.5304
## 1938 848.2953 768.9530 927.6376 848.2953 768.9530
## 1939 824.9833 745.6410 904.3256 824.9833 745.6410
## 1940 806.9241 727.5818 886.2664 806.9241 727.5818
## 1941 801.6043 722.2620 880.9466 801.6043 722.2620
## 1942 811.1336 731.7913 890.4759 811.1336 731.7913
## 1943 817.2703 737.9280 896.6126 817.2703 737.9280
## 1944 823.9198 744.5775 903.2621 823.9198 744.5775
## 1945 838.5404 759.1981 917.8827 838.5404 759.1981
## 1946 856.8138 777.4715 936.1561 856.8138 777.4715
## 1947 857.2625 777.9202 936.6048 857.2625 777.9202
## 1948 857.4448 778.1024 936.7871 857.4448 778.1024
## 1949 856.0162 776.6738 935.3585 856.0162 776.6738
## 1950 855.3676 776.0252 934.7100 855.3676 776.0252
## 1951 851.3491 772.0066 930.6916 851.3491 772.0066
## 1952 857.7762 778.4334 937.1189 857.7762 778.4334
## 1953 874.7876 795.4445 954.1306 874.7876 795.4445
## 1954 895.3786 816.0349 974.7223 895.3786 816.0349
## 1955 900.9243 821.5794 980.2692 900.9243 821.5794
## 1956 904.8085 825.4613 984.1556 904.8085 825.4613
## 1957 900.7923 821.4410 980.1437 900.7923 821.4410
## 1958 906.8756 827.5166 986.2347 906.8756 827.5166
## 1959 911.3900 832.0164 990.7636 911.3900 832.0164
## 1960 909.7148 830.3143 989.1153 909.7148 830.3143
## 1961 917.2558 837.8052 996.7064 917.2558 837.8052
## 1962 914.7993 835.2556 994.3430 914.7993 835.2556
## 1963 913.1991 833.4822 992.9159 913.1991 833.4822
## 1964 912.7858 832.7477 992.8239 912.7858 832.7477
## 1965 887.3445 806.7118 967.9773 887.3445 806.7118
## 1966 859.5042 777.7761 941.2323 859.5042 777.7761
## 1967 842.7083 758.9793 926.4372 842.7083 758.9793
## 1968 818.4888 731.1574 905.8202 818.4888 731.1574
## 1969 804.0474 710.3783 897.7165 804.0474 710.3783
## 1970 798.3679 693.9207 902.8152 798.3679 693.9207
## y.KS.CI.upr
## 1871 1216.1159
## 1872 1204.5270
## 1873 1192.5969
## 1874 1197.2451
## 1875 1194.1066
## 1876 1187.2398
## 1877 1175.6783
## 1878 1191.8924
## 1879 1196.7897
## 1880 1177.1730
## 1881 1153.4855
## 1882 1137.5165
## 1883 1133.5423
## 1884 1124.1436
## 1885 1119.6908
## 1886 1117.2195
## 1887 1122.3267
## 1888 1114.1022
## 1889 1128.8183
## 1890 1152.4355
## 1891 1169.5425
## 1892 1185.6959
## 1893 1191.7642
## 1894 1194.1759
## 1895 1183.4354
## 1896 1157.5245
## 1897 1117.8141
## 1898 1078.9281
## 1899 1030.2713
## 1900 998.8305
## 1901 975.1242
## 1902 953.5373
## 1903 949.4842
## 1904 938.6337
## 1905 930.3413
## 1906 936.6446
## 1907 937.2363
## 1908 953.9701
## 1909 956.5587
## 1910 942.3346
## 1911 917.7956
## 1912 893.9818
## 1913 878.7930
## 1914 897.0234
## 1915 914.6390
## 1916 945.2249
## 1917 951.0841
## 1918 934.7328
## 1919 920.6575
## 1920 914.1053
## 1921 908.8923
## 1922 909.6684
## 1923 909.0166
## 1924 905.0249
## 1925 897.4992
## 1926 901.6654
## 1927 903.6250
## 1928 913.3964
## 1929 926.8707
## 1930 921.6167
## 1931 924.4656
## 1932 933.5540
## 1933 941.5927
## 1934 951.3099
## 1935 954.0180
## 1936 946.0884
## 1937 935.2150
## 1938 927.6376
## 1939 904.3256
## 1940 886.2664
## 1941 880.9466
## 1942 890.4759
## 1943 896.6126
## 1944 903.2621
## 1945 917.8827
## 1946 936.1561
## 1947 936.6048
## 1948 936.7871
## 1949 935.3585
## 1950 934.7100
## 1951 930.6916
## 1952 937.1189
## 1953 954.1306
## 1954 974.7223
## 1955 980.2692
## 1956 984.1556
## 1957 980.1437
## 1958 986.2347
## 1959 990.7636
## 1960 989.1153
## 1961 996.7064
## 1962 994.3430
## 1963 992.9159
## 1964 992.8239
## 1965 967.9773
## 1966 941.2323
## 1967 926.4372
## 1968 905.8202
## 1969 897.7165
## 1970 902.8152
par(mfrow=c(2,2), cex=0.6)
# plot filtered state
plot.ts( cbind(y, y.KF.CI, Nile), plot.type="single",
col=c(1,4,4,4,1), lwd=c(2,2,2,2,1), lty=c(1,1,2,2,3), xlab="", ylab="", main="")
abline(v=1898)
legend("topright", legend = c("data","filtered state","90% confidence interval"),
col = c("black","blue","blue"), lty = c(1,1,2), lwd = c(1,1,1),bty="n", cex=0.9 )
plot( ts( c(y.out$v[-1]), start=1871), col="blue", lwd=2, xlab="", ylab="", main="forecast error")
abline(h=0, lty=3)
plot( ts( c(y.out$P)[-1], start=1871), col="blue", lwd=2, xlab="", ylab="", main="variance of state")
plot( ts( c(y.out$F)[-1], start=1871), col="blue", lwd=2, xlab="", ylab="", main="variance of forecast error")
par(mfrow=c(1,2), mar=c(3,3,2,1), cex=0.8)
# filtered
plot.ts( cbind(y, y.KF.CI, Nile), plot.type="single",
col=c(1,4,4,4,1), lwd=c(2,2,2,2,1), lty=c(1,1,2,2,3), xlab="", ylab="", main="")
abline(v=1898)
legend("topright", legend=c("data","filtered","90% confidence interval"), col=c(1,4,4), lty=c(1,1,2), lwd=2, bty="n", cex=0.9 )
# smoothed
plot.ts( cbind(y, y.KS.CI, Nile), plot.type="single",
col=c(1,2,2,2,1), lwd=c(2,2,2,2,1), lty=c(1,1,2,2,3), xlab="", ylab="", main="")
abline(v=1898)
legend("topright", legend=c("data","smoothed","90% confidence interval"), col=c(1,2,2), lty=c(1,1,2), lwd=2, bty="n", cex=0.9 )
