1 00:00:11,060 --> 00:00:16,390 In this lecture, we are going to apply auto arima to the airline passengers data center. 2 00:00:17,000 --> 00:00:21,860 Our main goal in this lecture will be to improve the results we got in the last lecture, which shouldn't 3 00:00:21,860 --> 00:00:24,080 be very difficult since we were just guessing. 4 00:00:24,560 --> 00:00:27,200 That is, if Auto Roraima works as it should. 5 00:00:29,150 --> 00:00:33,140 OK, so let's start by downloading the airline passengers data set once again. 6 00:00:37,360 --> 00:00:40,450 Next, let's import Pendas Nambiar matplotlib. 7 00:00:44,130 --> 00:00:50,160 Next, we're going to install the PMed Rhema Library, this is necessary because this library does not 8 00:00:50,160 --> 00:00:52,440 come preinstalled in Google CoLab. 9 00:00:59,790 --> 00:01:03,870 Next, we're going to load in our data set using that reads ESV. 10 00:01:06,410 --> 00:01:10,700 Next, we do a DFG head, as usual, to ensure our data looks as it should. 11 00:01:14,340 --> 00:01:19,590 Next, we're going to take the log of the passengers column and store this in the log passenger's column. 12 00:01:22,230 --> 00:01:27,200 This should be interesting to compare, since previously we found that logging the passengers column 13 00:01:27,390 --> 00:01:31,050 helps the model achieve better accuracy with the same Arima order. 14 00:01:32,980 --> 00:01:37,930 Next, we're going to split our data into train and test with the last 12 data points assigned to be 15 00:01:37,930 --> 00:01:38,330 the test. 16 00:01:42,830 --> 00:01:45,830 Next, we're going to import the PMed, a library. 17 00:01:49,660 --> 00:01:55,060 Next, we're going to call the auto Arima function, OK, so what are these arguments? 18 00:01:55,480 --> 00:01:58,400 The first argument is the data that one should be obvious. 19 00:01:58,900 --> 00:02:04,570 The second argument is, Trece, we're going to set this to true so that we can see which models auto 20 00:02:04,570 --> 00:02:07,610 Arima tests out while it tries to find the best model. 21 00:02:08,590 --> 00:02:14,200 The third argument is suppress warnings, which we are also going to set the true, as you know, stats 22 00:02:14,200 --> 00:02:15,400 models has lots of warnings. 23 00:02:15,400 --> 00:02:18,420 So this will make sure that those warnings are not printed out. 24 00:02:20,390 --> 00:02:26,720 Finally, we said seasonal equal to true, since we want to fit a full seasonal Arima model and we set 25 00:02:26,720 --> 00:02:30,410 em equal to 12 since we know that the seasonal period is 12. 26 00:02:31,220 --> 00:02:32,570 OK, so let's run this. 27 00:02:39,090 --> 00:02:43,020 All right, so once this is done, we call the function model that summary. 28 00:02:48,760 --> 00:02:52,660 So if you've ever worked with are, you should be very comfortable with this kind of output. 29 00:02:53,350 --> 00:02:53,740 All right. 30 00:02:53,740 --> 00:02:58,990 So basically, this is going to tell you the orders of the best model and other information, such as 31 00:02:58,990 --> 00:03:04,790 the AIC, the Bisi log likelihood and the P values of various statistical tests. 32 00:03:05,860 --> 00:03:10,600 You can also see the value of each parameter in the model and whether or not there are statistically 33 00:03:10,600 --> 00:03:11,320 significant. 34 00:03:12,330 --> 00:03:18,210 OK, so what's our best model, it looks like our best model is an R three for the non seasonal par 35 00:03:18,420 --> 00:03:21,550 and for the seasonal part, all it does is the first difference. 36 00:03:22,020 --> 00:03:24,590 Remember that this is the first seasonal difference. 37 00:03:24,750 --> 00:03:30,510 So instead of subtracting from one time step behind, it's subtracting from 12 times steps behind. 38 00:03:31,350 --> 00:03:35,430 You might actually want to try that in code to see if the result looks stationary. 39 00:03:36,430 --> 00:03:41,890 Also note that this is a much more parsimonious model than the non seasonal Arima we found earlier, 40 00:03:42,100 --> 00:03:44,170 which was in a Arima 12 one zero. 41 00:03:44,650 --> 00:03:47,830 This model has only three dates and one intercept. 42 00:03:53,880 --> 00:03:56,490 OK, so next, let's call model that predicts. 43 00:03:57,780 --> 00:04:03,570 As input will say, experience periods equals and test and return confidence equal to true so that we 44 00:04:03,570 --> 00:04:05,070 can plot the confidence interval. 45 00:04:09,320 --> 00:04:14,660 Next, we're going to plot our forecast along with the true values, so first we're going to create 46 00:04:14,660 --> 00:04:17,180 an axis object and set the figure size. 47 00:04:17,810 --> 00:04:20,910 Next, we're going to call the plot function on the test data frame. 48 00:04:21,170 --> 00:04:22,670 We'll give this low label data. 49 00:04:23,660 --> 00:04:26,890 Next, we're going to call the plot function on the test predictions. 50 00:04:27,290 --> 00:04:30,500 We'll give this the label forecast next. 51 00:04:30,510 --> 00:04:34,690 We're going to call the function of field between which will plot the confidence bounce. 52 00:04:35,240 --> 00:04:39,560 As with stats models, the confidence bounce are returned as a T by two array. 53 00:04:40,040 --> 00:04:46,070 So which row represents a timestamp in each column represents either the lower bound or the upper bound. 54 00:04:46,670 --> 00:04:49,810 Finally, we call X legend to plot the labels. 55 00:04:50,780 --> 00:04:52,160 OK, so let's run this. 56 00:04:56,170 --> 00:04:58,900 All right, and we see that our forecast is pretty good. 57 00:05:02,710 --> 00:05:07,700 Next, I want to plot a plot of the full model prediction, just like we did before. 58 00:05:08,260 --> 00:05:12,110 This will include both the forecast and the sample predictions. 59 00:05:12,640 --> 00:05:16,510 So first, we're going to need to call the function model, predict in sample. 60 00:05:16,960 --> 00:05:22,390 We'll pass in start equals zero and end equals minus one to get the full training predictions. 61 00:05:23,110 --> 00:05:28,990 Note that unlike stats models, we do not need to specify that we want the predictions to be accumulated. 62 00:05:29,200 --> 00:05:34,430 If you have differences in your model PMed Yarema, we'll just assume that automatically. 63 00:05:35,560 --> 00:05:38,750 Next, we call all the usual plot functions to make our plot. 64 00:05:39,700 --> 00:05:40,780 So let's run this. 65 00:05:45,480 --> 00:05:51,030 OK, so this pretty much looks as expected, both the training predictions and the forecasts look pretty 66 00:05:51,030 --> 00:05:51,360 good. 67 00:05:52,230 --> 00:05:58,260 One weird part about this graph is the beginning, where all we see is a flat line since we now have 68 00:05:58,260 --> 00:05:59,190 a seasonal model. 69 00:05:59,220 --> 00:06:04,980 It's not possible to make predictions for the first few times that this is due to the fact that we require 70 00:06:04,980 --> 00:06:08,220 values from a one season behind which don't exist. 71 00:06:12,930 --> 00:06:18,330 OK, so next, since we know that our model might perform better on the log data set, we're going to 72 00:06:18,330 --> 00:06:24,420 run our auto arima again, but on the log passengers, otherwise, all the arguments to this function 73 00:06:24,420 --> 00:06:25,650 are the same as before. 74 00:06:26,130 --> 00:06:27,270 So let's run this. 75 00:06:32,710 --> 00:06:35,160 All right, so now let's check model that summary. 76 00:06:42,550 --> 00:06:47,650 So as you can see, this finds a slightly different model, then the Ottawa rhema on the non logged 77 00:06:47,650 --> 00:06:48,260 data set. 78 00:06:48,910 --> 00:06:53,860 Now we have two auto regressive coefficients and one seasonal moving average coefficient. 79 00:06:54,430 --> 00:06:58,240 Again, we do seasonal differences rather than normal seasonal difference in. 80 00:07:03,520 --> 00:07:09,220 All right, so again, we have a plot of the forecast against the true values, I'm not going to explain 81 00:07:09,220 --> 00:07:12,100 this code since it's exactly the same as the previous code. 82 00:07:18,200 --> 00:07:21,200 So qualitatively, the results still look pretty good. 83 00:07:25,370 --> 00:07:30,950 Next, let's look at the full train and test predictions again, this is the same code as before. 84 00:07:36,120 --> 00:07:41,220 OK, so the plot is kind of hard to see, since the first few values stretch all the way down to zero, 85 00:07:41,920 --> 00:07:44,850 you could block this if you wanted, but I'm not going to bother. 86 00:07:49,210 --> 00:07:54,400 Now, since both of these models look to be performing pretty well, it might be a good idea to test 87 00:07:54,400 --> 00:07:57,900 them quantitatively, that is to use the IMC. 88 00:07:58,570 --> 00:08:01,860 You've seen this function before, so there's no need to explain it again. 89 00:08:05,970 --> 00:08:07,170 So let's run this. 90 00:08:08,960 --> 00:08:14,600 All right, so this is kind of the opposite result from what we saw previously in this case, the non 91 00:08:14,810 --> 00:08:16,190 diversion performs better. 92 00:08:17,060 --> 00:08:23,360 Interestingly, we get an army of about 18, something which is on par with our previous Arima 12 one 93 00:08:23,360 --> 00:08:26,030 zero model, the one we found by just guessing. 94 00:08:30,050 --> 00:08:34,530 Now, the next thing I would like to do in the script is to test Otto Arima itself. 95 00:08:35,090 --> 00:08:39,500 I want to know, can Otto Arima find the model that we found last time? 96 00:08:40,010 --> 00:08:43,900 So I'm going to call Otto again on the log passengers. 97 00:08:44,510 --> 00:08:48,010 This time I'm going to set a few more parameters into this function. 98 00:08:48,740 --> 00:08:55,520 First, I'm gonna set Max P equal to twelve and Max Q equal to two if you look at the documentation. 99 00:08:55,590 --> 00:08:57,950 These values are set quite low by default. 100 00:08:58,310 --> 00:09:03,530 So if we want auto arima to even have a chance of finding the model that we found, we have to set them 101 00:09:03,530 --> 00:09:04,130 higher. 102 00:09:05,210 --> 00:09:08,230 Note that we also have to set a parameter called Max Order. 103 00:09:08,840 --> 00:09:15,120 This value has a limit on the total sum of all the Arima orders, which promotes finding a simpler model. 104 00:09:15,890 --> 00:09:18,900 Again, the default value is five, which is very small. 105 00:09:19,160 --> 00:09:23,470 So in order to have a chance of finding our previous model, we need to set this higher. 106 00:09:24,050 --> 00:09:29,150 I've set it to 14 since that's the sum of 12 into next. 107 00:09:29,150 --> 00:09:31,160 We have an argument called stepwise. 108 00:09:31,760 --> 00:09:37,160 This controls whether we use an intelligent algorithm to search the hyper parameter space or if you 109 00:09:37,160 --> 00:09:41,250 use a naive grid, search the default value for this is true. 110 00:09:41,480 --> 00:09:44,570 So for now, we'll leave it like this and see how it performs. 111 00:09:45,290 --> 00:09:50,750 Finally, we set Seasonale equal to false, since our goal right now is to compare these results to 112 00:09:50,750 --> 00:09:52,430 the results from the previous script. 113 00:09:53,330 --> 00:09:54,770 OK, so let's run this. 114 00:10:00,920 --> 00:10:05,560 OK, so what's interesting about this is that it doesn't even try that many values. 115 00:10:05,840 --> 00:10:12,120 In fact, we never get past P equals five, but perhaps this is a better model than what we found before. 116 00:10:19,920 --> 00:10:25,770 If we do a model that summary, we can see that the presumably best model is a rhema for one two. 117 00:10:30,840 --> 00:10:33,000 Next, let's plot our model forecast. 118 00:10:37,820 --> 00:10:40,430 OK, so we can see that this forecast is pretty bad. 119 00:10:44,690 --> 00:10:47,840 Next, let's plot both the terrain and test predictions together. 120 00:10:52,760 --> 00:10:54,010 Again, not so great. 121 00:10:58,030 --> 00:10:59,890 Finally, let's check the arms, see? 122 00:11:02,850 --> 00:11:06,540 So as you can see, it's pretty high relative to what we've seen before. 123 00:11:07,320 --> 00:11:13,220 OK, so at this point what I want to do is to compare what we just did to a full grid search. 124 00:11:13,800 --> 00:11:14,910 So let's go back up. 125 00:11:21,340 --> 00:11:23,500 And let's set stepwise false. 126 00:11:27,510 --> 00:11:28,950 And let's run this all again. 127 00:11:33,560 --> 00:11:39,350 OK, so this time we can see that auto arima is actually testing out all the possible combinations of 128 00:11:39,350 --> 00:11:46,580 models, not that it would not be a good idea to do this if you had a full screamo model that is both 129 00:11:46,580 --> 00:11:48,830 the Nazis in a part and the seasonal part. 130 00:11:49,460 --> 00:11:52,940 Again, this goes back to the concept of the curse of dimensionality. 131 00:11:53,660 --> 00:11:57,890 Currently, we're only testing out two hyper parameters, which is a two dimensional grid. 132 00:11:58,520 --> 00:12:03,710 If you have on average 10 values of each type of parameter, then that would be one hundred values to 133 00:12:03,710 --> 00:12:04,340 test. 134 00:12:04,940 --> 00:12:09,800 If you had the seasonal part to, then you would have 10 to the power for values to test. 135 00:12:10,100 --> 00:12:13,490 That's 10000, which is 100 times more than one hundred. 136 00:12:18,280 --> 00:12:21,050 OK, so in any case, let's run model that summary. 137 00:12:23,350 --> 00:12:27,410 As you can see, the best model appears to be a rhema 12 one one. 138 00:12:27,880 --> 00:12:32,140 This is very close to our manually found model, which was 12 one zero. 139 00:12:36,840 --> 00:12:39,120 OK, so next, let's plot the forecast. 140 00:12:43,220 --> 00:12:44,640 OK, so that looks pretty good. 141 00:12:46,070 --> 00:12:51,890 So already we can see that perhaps auto Arima might not necessarily find the best model with the default 142 00:12:51,890 --> 00:12:52,690 settings. 143 00:12:57,600 --> 00:13:00,530 Next, let's plot the terrain and test predictions together. 144 00:13:03,690 --> 00:13:05,280 Again, this looks pretty good. 145 00:13:09,770 --> 00:13:12,080 Finally, let's print out the pharmacy. 146 00:13:15,260 --> 00:13:17,400 OK, so we get about twenty three something. 147 00:13:18,110 --> 00:13:23,200 This is actually quite interesting because this is a worst model than the model we found on our own. 148 00:13:23,780 --> 00:13:28,490 As you recall, a remote 12 one zero achieves 18 something or AMZI. 149 00:13:29,210 --> 00:13:32,210 Furthermore, that model has less parameters than this one. 150 00:13:32,570 --> 00:13:36,530 So this model actually gets a higher penalty for having one more parameter. 151 00:13:37,250 --> 00:13:43,790 So what this tells us is that finding the best AIC is not the same as minimizing the out of sample error. 152 00:13:44,630 --> 00:13:50,570 They are similar ideas in concept, but one is looking for a simpler model purely to look for a simpler 153 00:13:50,570 --> 00:13:55,580 model, while the other doesn't actually care about model complexity, but rather accuracy. 154 00:13:56,150 --> 00:13:59,150 So these two ideas are related, but they are not the same.