Lo scopo di questo progetto è quello di analizzare diversi algoritmi per le previsioni finanziarie di alcune azioni molto famose. Questa analisi è stata condotta da due aspiranti Data Scientist, pertanto non costituisce e non è da considerarsi una consulenza finanziaria.

Esistono molti modelli statistici per i mercati finanziari. Il comportamento di questi modelli è simile a quello dei modelli dei fenomeni naturali. Cioè, entrambi sono influenzati da variabili sconosciute e instabili. Ciò comporta un’elevata e imprevedibile volatilità. È quindi quasi impossibile prevedere il comportamento futuro.

** Per l’analisi delle serie storiche abbiamo consultato questo PDF del Prof. Vito Ricci: –> https://cran.r-project.org/doc/contrib/Ricci-ts-italian.pdf **

Per il trattamento delle classi (ts, xts, zoo) abbiamo consultato questo PDF di Rotman School University of Toronto –> https://tdmdal.github.io/r-tutorial-202021-winter/r_timeseries_finance_pkgs.pdf

Librerie usate:

require(xts)
require(doParallel)
require(caret)
require(dynlm)
library(zoo)
library(nnet)
library(MASS)
library(Rsolnp)
library(nlme)
library(TTR)
library(depmixS4)
library(ggplot2)
library(reshape2)
library(PerformanceAnalytics)
library(quantmod)
library(forecast)
library(highcharter)
library(tseries)
library(timeSeries)
library(TSstudio)
library(dplyr)
library(ggthemes)
require(kableExtra)
library(dynlm)
library(tidyr)
library(TSA)
library(gridExtra)
library(corrplot)
library(plotly)
library(GGally)
library(data.table)
library(lubridate)
library(DT)
library(e1071)
library(doParallel)
library(knitr)
library(ggdark)
library(ggplot2)
library(depmixS4)
# library(tensorflow)
# library(keras)
# options(kableExtra.latex.load_packages = FALSE)
# ##### NOTA BENE (SE HAI UN MAC CON CPU SILICON GLI ULTIMI AGIONAMENTI DI TENSOR E KERAS NON FUNZIONANO)--> Ho dovuto eseguire il downgrade di tensorflow-macos alla versione 2.9.0 e tensorflow-metal alla 0.5.0 per eliminare tali errori. Altri su questo sito hanno segnalato problemi simili con versioni più recenti. Apple deve ancora fornire versioni funzionanti per le nuove versioni di TF.
# #https://developer.apple.com/forums/thread/722361

Utilizziamo la libreria “quantmod” per scaricare gli indici azionari interessati

price_apple <- getSymbols("AAPL", auto.assign=FALSE, from="2018-01-01", to="2022-12-31")
price_google <- getSymbols("GOOGL",auto.assign=FALSE,from="2018-01-01",to="2022-12-31")
price_amazon <- getSymbols("AMZN",auto.assign=FALSE,from="2018-01-01",to="2022-12-31")
price_tesla <- getSymbols("TSLA",auto.assign=FALSE,from="2018-01-01",to="2022-12-31")

Con la libreria “highcharter” visualizziamo gli indici azionari interessati:

highchart(type="stock") %>% 
  hc_add_series(Cl(price_apple), name="AAPL") %>% 
  hc_add_series(Cl(price_google), name="GOOGL") %>% 
  hc_add_series(Cl(price_amazon), name="AMZN") %>% 
  hc_add_series(Cl(price_tesla), name="TSLA") %>% 
  hc_title(text="<b>AAPL vs GOOGL vs AMZN Closing Price</b>")

Da questo primo grafico vediamo che il prezzo più alto c’è l’ha Apple

Calcoliamo i rendimenti azionari

return_apple <- dailyReturn(Cl(price_apple))
return_google <- dailyReturn(Cl(price_google))
return_amazon <- dailyReturn(Cl(price_amazon))
return_tesla <- dailyReturn(Cl(price_tesla))

Creiamo una nuova variabile (returns) che conterrà tutti i rendimenti azionari

returns <- data.frame(return_apple,return_google,return_amazon,return_tesla)
names(returns) <- c("return_apple","return_google","return_amazon","return_tesla")
returns <- as.xts(returns)
charts.PerformanceSummary(returns,main="Rendimenti giornalieri Apple vs Google vs Amazon vs Tesla 2018-2022")

TSLA presenta un un trend Ribassista nei mesi di Aprile e di settembre 2022 e APPL presenta un andamento rialzista a medio lungo termine, tuttavia si sono rivelati avere il più alto rendimento cumulativo rispetto a gli altri. Nella parte di drawdown possiamo vedere che ad oggi APPL è posizionata al disopra di tutte le altre. Questo significa che ha avuto un calo solo a marzo 2020 (come tutte le altre) e fin da allora non è mai più scesa a differenza delle altre

Nella parte di drawdown possiamo vedere che ad oggi APPL è posizionata al disopra di tutte le altre. Questo significa che ha avuto un calo solo a marzo 2020 (come tutte le altre) e fin da allora non è mai più scesa a differenza delle altre.

se vuoi sapere di più sul drawdown puoi cliccare qui –>https://www.avatrade.it/education/market-terms/what-is-drawdown

Dopo aver fatto una breve analisi, sugli indici azionari scelti, abbiamo deciso di continuare ad approfondire la nostra analisi basandoci solo l’indice AAPL

chartSeries(price_apple, theme="black",
            TA="addVo();addBBands();addCCI()", subset = '2018-01::')

** Questo grafico ci mostra l’andamento dell’indice apple con le bande di bollinger, il volume e il CCI (Commodity Channel Index) se vuoi ulteriori informazioni sulle bande di bollinger e il CCI puoi visitare questi link: https://www.avatrade.it/education/technical-analysis-indicators-strategies/cci-trading-strategies https://www.ig.com/it/strategie-di-trading/cosa-sono-le-bande-di-bollinger-e-come-usarle-nel-trading-190122**

highchart(type="stock") %>% 
  hc_add_series(price_apple) %>% 
  hc_add_series(SMA(na.omit(Cl(price_apple)),n=50),name="SMA(50)") %>% 
  hc_add_series(SMA(na.omit(Cl(price_apple)),n=200),name="SMA(200)") %>% 
  hc_title(text="<b>APPLE Price Candle Stick Chart 2018-2022</b>") %>% 
  hc_add_theme(hc_theme_darkunica())

Tali visualizzazioni possono essere molto utili quando è necessario vedere più dettagliatamente i dati. Le Medie mobili singole (SMA) sono anche più facili da vedere con un grafico interattivo. rispetto ai grafici statici. Questo grafico può essere utilizzato per evidenziare importanti analisi tecniche che influenzano il nostro processo decisionale. come la “croce d’oro” o la “croce della morte”. Se vuoi saperne di più su quest’ultima puoi dare un occhiata a questo sito > https://www.mazzieroresearch.com/golden-cross-che-cose/

“ts_decompose” è una funzione della libreria TS_studio, ci permette di vedere quattro grafici interattivi sui prezzi di chiusura

closing_pr <- Cl(to.monthly(price_apple))
dc_ts <- ts_decompose(as.ts(closing_pr, start=c(2018)))
dc_ts
dc <- decompose(as.ts(closing_pr, start=c(2018,1,1)))
dc$seasonal
##             Jan        Feb        Mar        Apr        May        Jun
## 2018  5.7966447 -3.3755472 -2.8190606 -3.4712485 -9.5817688 -7.1869758
## 2019  5.7966447 -3.3755472 -2.8190606 -3.4712485 -9.5817688 -7.1869758
## 2020  5.7966447 -3.3755472 -2.8190606 -3.4712485 -9.5817688 -7.1869758
## 2021  5.7966447 -3.3755472 -2.8190606 -3.4712485 -9.5817688 -7.1869758
## 2022  5.7966447 -3.3755472 -2.8190606 -3.4712485 -9.5817688 -7.1869758
##             Jul        Aug        Sep        Oct        Nov        Dec
## 2018  1.8963295  8.5080738  0.8137514 -0.3231225  2.4619044  7.2810195
## 2019  1.8963295  8.5080738  0.8137514 -0.3231225  2.4619044  7.2810195
## 2020  1.8963295  8.5080738  0.8137514 -0.3231225  2.4619044  7.2810195
## 2021  1.8963295  8.5080738  0.8137514 -0.3231225  2.4619044  7.2810195
## 2022  1.8963295  8.5080738  0.8137514 -0.3231225  2.4619044  7.2810195

L’output ci mostra quattro grafici dei nostri dati sui prezzi di chiusura, che sono:

Observeded <- ci spiega l’andamento dell’indice azionario

Trend <- con il trend possiamo vedere la significativa tendenza al rialzo, iniziata intorno alla metà dell’anno 2019.

**Seasonal <- fluttuazione stagionale ripetitiva dei dati. Il prezzo di chiusura di APPL tende a raggiungere il più alto ad Agosto (Generalmente nel mese di settembre vengono presentati nuovi dispositivi come iphone e computer, questo potrebbe spiegare l’andamento rialzista) e il più basso a Maggio. Guardando questo modello, possiamo dire che nel complesso, il momento giusto per vendere questo titolo era nel mese di Dicembre, Gennaio ma soprattutto Agosto e il momento giusto per acquistare era nel mese di Maggio e Giugno, questo lo possiamo vedere meglio eseguendo questa riga di codice —> dc$seasonal.**

**Random <- fluttuazione irregolare o casuale non catturata dal trend e dalla stagionalità, L’attuale situazione della pandemia di Covid-19 è un esempio del fatto che potrebbe causare questa fluttuazione casuale.**

seas <- ts(price_apple, start=c(2018), end=c(2022, 12), frequency=12)

seas
##          AAPL.Open AAPL.High AAPL.Low AAPL.Close AAPL.Volume AAPL.Adjusted
## Jan 2018   42.5400   43.0750  42.3150    43.0650   102223600      40.95049
## Feb 2018   43.1325   43.6375  42.9900    43.0575   118071600      40.94336
## Mar 2018   43.1350   43.3675  43.0200    43.2575    89738400      41.13354
## Apr 2018   43.3600   43.8425  43.2625    43.7500    94640000      41.60185
## May 2018   43.5875   43.9025  43.4825    43.5875    82271200      41.44733
## Jun 2018   43.6375   43.7650  43.3525    43.5825    86336000      41.44259
## Jul 2018   43.2900   43.5750  43.2500    43.5725    95839600      41.43308
## Aug 2018   43.6475   43.8725  43.6225    43.8200    74670800      41.66843
## Sep 2018   44.0450   44.3400  43.9125    44.2725   101672400      42.09871
## Oct 2018   44.4750   44.8475  44.0350    44.0475   118263600      41.88477
## Nov 2018   44.0375   44.8125  43.7675    44.7750   137547200      42.57653
## Dec 2018   44.8425   45.0250  44.5625    44.8150   124773600      42.61456
## Jan 2019   44.6525   44.8950  44.3525    44.6150   129700400      42.42439
## Feb 2019   44.3250   44.4450  44.1500    44.2500   108434400      42.07732
## Mar 2019   44.3250   44.8600  44.2050    44.2600   130756400      42.08682
## Apr 2019   44.3125   44.3250  43.3000    43.5550   204420400      41.41644
## May 2019   43.6275   43.7375  42.6325    42.7775   166116000      40.67711
## Jun 2019   43.0000   43.0000  42.5150    42.8775   156572000      40.77220
## Jul 2019   42.5400   42.5400  41.7675    41.9900   202561600      39.92828
## Aug 2019   41.3825   41.8425  41.1750    41.7425   184192800      39.69294
## Sep 2019   41.7175   42.1100  41.6250    41.8575   129915600      39.80228
## Oct 2019   41.7925   42.1550  41.6900    41.9450   188923200      39.88549
## Nov 2019   41.5000   41.7000  40.0250    40.1250   346375200      38.15486
## Dec 2019   39.7750   40.9700  39.0000    39.1225   290954000      37.20157
## Jan 2020   38.7075   40.9300  38.5000    40.7575   272975200      38.75629
## Feb 2020   40.7725   40.8500  39.7675    39.8850   206434400      37.92663
## Mar 2020   40.0725   40.2500  38.7575    38.7875   217562000      36.88302
## Apr 2020   39.2675   39.4725  37.5600    39.1025   282690400      37.33416
## May 2020   39.6250   40.9725  39.3775    40.6775   243278000      38.83792
## Jun 2020   40.4875   41.1875  40.4125    41.0850   130196800      39.22700
## Jul 2020   40.7600   41.8850  40.7200    41.8425   162579600      39.95025
## Aug 2020   42.4475   43.2725  42.2500    43.2475   204588800      41.29171
## Sep 2020   43.0900   43.7050  42.9425    43.1075   160704400      41.15802
## Oct 2020   43.0125   43.5650  42.8550    42.9625   135722000      41.01959
## Nov 2020   43.2075   43.5300  42.7525    42.7675   149886400      40.83342
## Dec 2020   42.9500   43.4875  42.9275    43.1250   123967600      41.17475
## Jan 2021   43.4175   43.9125  43.3850    43.8750   135249600      41.89082
## Feb 2021   44.0875   44.8475  44.0525    44.7425   152648800      42.71910
## Mar 2021   44.7750   45.1200  44.5400    44.5975   155712400      42.58064
## Apr 2021   44.8150   45.1550  44.5125    44.5300   151128400      42.51620
## May 2021   44.6350   44.9450  43.1650    43.7500   195208000      41.77147
## Jun 2021   43.2000   44.0750  43.1125    44.0525   153816000      42.06031
## Jul 2021   43.8025   44.4350  43.6300    44.2050   113605600      42.20590
## Aug 2021   44.4775   44.5625  44.0325    44.1675    95154000      42.17009
## Sep 2021   43.7350   43.9625  43.5675    43.7575   126814000      41.77864
## Oct 2021   43.8700   44.2800  43.7675    44.2350    95096400      42.23455
## Nov 2021   44.4900   45.0000  44.3475    44.9950   128740800      42.96017
## Dec 2021   45.0725   45.5975  45.0525    45.4300   128828400      43.37550
## Jan 2022   45.6475   45.8750  44.8100    44.9925   126774000      42.95779
## Feb 2022   45.0800   45.1300  44.4525    44.6100   117473600      42.59260
## Mar 2022   44.6250   45.0600  44.5175    44.6625    90975200      42.64271
## Apr 2022   44.6625   44.7800  44.4050    44.5050   157618800      42.49233
## May 2022   44.3300   44.3675  43.4150    43.8250   133787200      41.84309
## Jun 2022   43.8100   44.2000  43.7350    43.8100    78597600      41.82877
## Jul 2022   43.7600   43.7725  42.8150    42.8175   148219600      40.88115
## Aug 2022   42.5000   43.1700  42.1500    42.2125   165963200      40.30351
## Sep 2022   42.0975   42.4800  41.2350    41.2350   164115200      39.37021
## Oct 2022   42.0175   43.2750  41.6100    43.1925   150164800      41.23919
## Nov 2022   43.4200   43.7875  41.7300    42.0850   163690400      40.18178
## Dec 2022   41.8125   42.5050  41.2975    41.6200   166674000      39.73780
## attr(,"src")
## [1] yahoo
## attr(,"updated")
## [1] 2023-01-08 15:00:42 CET
## attr(,"index")
##  [1] 1514851200 1514937600 1515024000 1515110400 1515369600 1515456000
##  [7] 1515542400 1515628800 1515715200 1516060800 1516147200 1516233600
## [13] 1516320000 1516579200 1516665600 1516752000 1516838400 1516924800
## [19] 1517184000 1517270400 1517356800 1517443200 1517529600 1517788800
## [25] 1517875200 1517961600 1518048000 1518134400 1518393600 1518480000
## [31] 1518566400 1518652800 1518739200 1519084800 1519171200 1519257600
## [37] 1519344000 1519603200 1519689600 1519776000 1519862400 1519948800
## [43] 1520208000 1520294400 1520380800 1520467200 1520553600 1520812800
## [49] 1520899200 1520985600 1521072000 1521158400 1521417600 1521504000
## [55] 1521590400 1521676800 1521763200 1522022400 1522108800 1522195200
## attr(,"index")attr(,"tzone")
## [1] UTC
## attr(,"index")attr(,"tclass")
## [1] Date
ts_seasonal(seas, type = "all") 
## Warning in ts_seasonal(seas, type = "all"): The input object is a 'mts' class,
## by defualt will use only the first series as an input

ts_seasonal è anche una funzione della libreria “TS_studio”, in questo caso mettendo type=all , possiamo vedere dal grafico la stagionalità per anni e per mesi, ci consente di identificare i modelli stagionali

Modelli markoviani **Il Markov nascosto è modellato da un insieme predeterminato di gaussiane. Il problema del rilevamento del dominio è un problema di apprendimento non supervisionato, poiché il numero di stati non è noto a priori e non esiste una “verità di base” per “addestrare” l’HMM.

Questa sezione risolve due problemi di modellazione distinti: il primo consiste nell’adattare un HMM con due stati di modalità ai rendimenti dell’indice AAPL Il secondo consiste nell’utilizzare tre stati; Abbiamo usato la libreria “DepmixS4”**

L’oggetto della serie temporale applerest può essere tracciato, mostrando i periodi volatili intorno al 2018 e al 2023

plot(appleRest)

Un modello di Markov nascosto a due stati viene montato utilizzando l’algoritmo EM “Expectation-maximization”

breve guida –> https://rstudio-pubs-static.s3.amazonaws.com/154174_78c021bc71ab42f8add0b2966938a3b8.html

#Adattiamo un modello di Markov nascosto con due stati:
hmm <- depmix(returns ~ 1, family = gaussian(), nstates = 2, data=data.frame(returns=returns))
hmmfit <- fit(hmm, verbose = FALSE)
## converged at iteration 52 with logLik: 3213.266
post_probs <- posterior(hmmfit)
## Warning in .local(object, ...): Argument 'type' not specified and will default
## to 'viterbi'. This default may change in future releases of depmixS4. Please
## see ?posterior for alternative options.
plot(returns, type='l', main='Regime Detection', xlab='', ylab='Returns')

matplot(post_probs[,-1], type='l', main='Regime Posterior Probabilities', ylab='Probability')
legend(x='bottomleft', c('Regime #1','Regime #2'), fill=1:2, bty='n')

Traccia il flusso di returns e il regime posteriore probabilità dei regimi separati In questo caso abbiamo un HMM con 2 stati nascosti, questo grafico ci mostra due curve di probabilità, una per ogni stato nascosto. L’asse x del grafico rappresenta il tempo e l’asse y rappresenta la probabilità di trovarci in un particolare stato nascosto. l’altezza della curva ad ogni passo temporale, indica la probabilità di essere in quello stato in quel momento e la larghezza la rappresenta la durata dello stato

HMM con 3 stati

plot(returns, type='l', main='Regime Detection', xlab='', ylab='Returns')

matplot(post_probs1[,-1], type='l', main='Regime Posterior Probabilities', ylab='Probability')
legend(x='bottomleft', c('Regime #1','Regime #2', 'Regime #3'), fill=1:3, bty='n')

In questo caso abbiamo 3 curve, come detto precedentemente ognuna delle quali rappresenta la probabilità di essere in uno dei 3 stati.

Navigando su Medium abbiamo Letto un articolo interessante che tratta Hidden Markov Models for Time Series in R studio [Stock Market Data] scritto da Ankit Sekseria, abbiamo provato a replicarlo, e adattare il suo modello ai nostri dati. **link* https://medium.com/analytics-vidhya/hidden-markov-models-for-time-series-in-r-studio-5ae2b9fb0701

Estraiamo informazioni relative all’indice AAPL, il set di dati contiene informazioni relative al prezzo di apertura al 01-01-2018 e chiusura al 02-01-2023 delle azioni, al prezzo alto e basso e al volume delle azioni. Siamo interessati a modellare con il nostro HMM la differenza tra il valore di chiusura e il valore di apertura del giorno corrente.

mod1 <- depmix(AAPL.Close ~ 1, family = gaussian(), nstates = 5,
               data = AAPL_train)
set.seed(1)
fm2 <- fit(mod1, verbose = FALSE) 
probs <- posterior(fm2)
## Warning in .local(object, ...): Argument 'type' not specified and will default
## to 'viterbi'. This default may change in future releases of depmixS4. Please
## see ?posterior for alternative options.
head(probs)
##   state         S1        S2           S3            S4           S5
## 1     2 0.00000000 1.0000000 0.000000e+00  0.000000e+00 0.000000e+00
## 2     2 0.03625842 0.9637416 0.000000e+00 1.304355e-110 3.309715e-55
## 3     2 0.03339146 0.9665772 4.250485e-58  3.136680e-05 1.804330e-38
## 4     2 0.03975022 0.9602045 3.038961e-41  3.398249e-05 1.129905e-05
## 5     2 0.03440663 0.9655479 1.485854e-08  3.680699e-05 8.606687e-06
## 6     2 0.03566405 0.9642929 1.169993e-08  3.355836e-05 9.447186e-06
AAPL_predict <- cbind(AAPL_subset$AAPL.Close, probs$state)

head(AAPL_predict)
##            AAPL.Close probs.state
## 2018-01-02    43.0650           2
## 2018-01-03    43.0575           2
## 2018-01-04    43.2575           2
## 2018-01-05    43.7500           2
## 2018-01-08    43.5875           2
## 2018-01-09    43.5825           2
chartSeries(AAPL_predict [,1])

addTA(AAPL_predict[AAPL_predict [,2]==1,1],on=1, type= "p", col=5,pch=25)

addTA(AAPL_predict[AAPL_predict [,2]==2,1],on=1, type= "p", col=6,pch=24)

addTA(AAPL_predict[AAPL_predict [,2]==3,1],on=1, type= "p", col=7,pch=23)

addTA(AAPL_predict[AAPL_predict [,2]==4,1],on=1, type= "p", col=8,pch=22)

addTA(AAPL_predict[AAPL_predict [,2]==5,1],on=1, type= "p", col=10,pch=21)

Arima Forecasting

** L’idea alla base delle previsioni delle serie temporali è quella di utilizzare i dati passati e attuali per prevedere il futuro. Anche se non tutti i modelli utilizzati per le previsioni possono essere accurati al 100%, i risultati sono comunque generalmente molto utili per prendere decisioni sul futuro. Esistono molti modelli e algoritmi per la previsione dei dati delle serie temporali.

Come analisti, dobbiamo decidere quale rapporto di separazione ci permetterà di presentare i dati di addestramento e di test senza uneccessivo costocomputazionale perl’addestramento del modello. Il rapporto più comune è 80% train - 20% test, ma dipende dai dati e dallo scopo dellaprevisione. *splitting:

*train data: dati utilizzati per adattarsi al modello

*test data: dati utilizzati per valutare il modello

In questo caso vogliamo prevedere i prezzi dei titoli su 100 giorni nel set di dati. Utilizzeremo quindi le ultime 100 osservazioni come dati di prova e il resto dei dati come dati di addestramento.***

# Numero di giorni che vogliamo prevedere
n <- 100

Effettuiamo lo split del dataset dividendo i dati in train e test:

# Splitting the data
train <- head(Cl(price_apple), length(Cl(price_apple))-n)
test <- tail(Cl(price_apple), n)

Il metodo ingenuo (Naive Bayes) è un metodo di previsione in cui l’ultima osservazione viene utilizzata come risultato previsto dei dati. Viene utilizzato come modello di base per le previsioni. Se i risultati del nostro modello sono peggiori di quelli del metodo ingenuo, decidiamo di non utilizzare questo modello.

fc_na <- naive(train, h=n)

plot dei risultati:

autoplot(fc_na) +
  autolayer(ts(test, start=length(train)), series = "forecast")

La linea blu è la media della nostra previsione, mentre le aree più scure e chiare più scure, che rappresentano rispettivamente gli intervalli di confidenza dell’80% e del 95%. Se confrontiamo il risultato con i dati di prova effettivi, possiamo vedere che ci sono differenze tra loro.

** Arima Model **

**Il modello autoregressivo integrato a media mobile (ARIMA) è una combinazione di un modello autoregressivo, un modello di integrazione per differenza e un modello a media mobile. Il modello autoregressivo (AR) determina la relazione tra un’osservazione e diverse osservazioni ritardate; per utilizzare il modello ARIMA i dati delle serie temporali devono essere stazionari, il che può essere ottenuto differenziando i dati. Infine, il modello di media mobile (MA) definisce la relazione tra l’errore residuo del modello di media mobile rispetto alle osservazioni ritardate e i valori osservati. I modelli ARIMA sono solitamente indicati come ARIMA(p,d,q), dove p, d e q sono parametri stimati positivamente; l’ordine AR (p) è il numero di osservazioni ritardate nel modello. Per d, d è il numero di volte in cui i dati effettivi vengono differenziati fino allo stato stazionario. In genere, non ci sono più di due differenziali prima di raggiungere lo stato stazionario; l’ordine MA (q) è la dimensione della finestra della media mobile. Esistono modi per determinare valori appropriati per questi parametri nel nostro modello, ma è difficile. Fortunatamente, R ha una funzione auto.arima() che fa questo lavoro per noi. Utilizzando questa funzione, possiamo ottenere due tipi di ARIMA: ARIMA non stagionale (che abbiamo appena descritto) e ARIMA stagionale. Ovviamente, l’ARIMA non stagionale non include la parte stagionale dei dati, mentre l’ARIMA stagionale sì. Vediamo quale fornisce i migliori risultati di previsione per i nostri dati.

**ARIMA non stagionale**
model_non <- auto.arima(train, seasonal=FALSE)
fc_non <- forecast(model_non, h=n)

plot dei risultati con Arima non stagionale

autoplot(fc_non)+
  autolayer(ts(test, start= length(train)), series="forecast")

I risultati ottenuti utilizzando l’ARIMA non stagionale hanno mostrato una tendenza al rialzo. Tuttavia, ci sono ancora alcune differenze tra i due quando si confrontano con i dati dei test reali.

** Arima stagionale**
model_s <- auto.arima(train)
fc_s <- forecast(model_s, h=n)
autoplot(fc_s)+
  autolayer(ts(test, start= length(train)), series="forecast")

I risultati mostrano che auto.arima() dà gli stessi risultati indipendentemente dal fatto che la parte stagionale dei dati sia inclusa o meno. Ciò significa che la parte stagionale dei dati non è significativa e che l’ARIMA non stagionale è il miglior modello ARIMA basato su auto.arima() per il nostro set di dati.

 **Forecast Evaluation - valutazione delle previsioni**

Dopo aver previsto i dati, l’ultima fase consiste nel valutare i risultati delle previsioni. La valutazione delle previsioni viene effettuata controllando che i residui siano coerenti con le ipotesi sui residui e confrontando gli indicatori di accuratezza. Controllo dei residui: Il residuo è la differenza tra i dati previsti e quelli effettivi. Un buon modello è quello in cui i residui sono distribuiti in modo casuale e non c’è un modello evidente. Questa sezione descrive le ipotesi residue che devono essere soddisfatte. La distribuzione è normale (media = 0) e può essere convalidata da una curva normale. La curva normale deve essere a campana. Ha una varianza costante, che può essere confermata da un grafico dei residui. La varianza costante è indicata da fluttuazioni costanti nei dati. Non c’è autocorrelazione, il che può essere confermato dal grafico ACF e dal test di Ljung-Box. L’autocorrelazione può essere rilevata sul grafico ACF se c’è una linea che supera il limite superiore o inferiore. Nel test di Ljung-Box il valore p deve essere maggiore di 0,05 per soddisfare questa ipotesi; se i risultati del grafico ACF e del test di Ljung-Box differiscono, si raccomanda di utilizzare i risultati del test di Ljung-Box.

residual of naive method

checkresiduals(fc_na)

## 
##  Ljung-Box test
## 
## data:  Residuals from Naive method
## Q* = 28.123, df = 10, p-value = 0.001725
## 
## Model df: 0.   Total lags used: 10

residual of arima model

checkresiduals(fc_non)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1) with drift
## Q* = 19.183, df = 9, p-value = 0.02368
## 
## Model df: 1.   Total lags used: 10

I risultati mostrano che il metodo ingenuo non soddisfa l’ipotesi di assenza di autocorrelazione, mentre l’ARIMA soddisfa tutte le ipotesi. Pertanto, sulla base del controllo dei residui, il nostro modello ARIMA ha fornito risultati migliori rispetto al metodo ingenuo.

Confronto tra gli indicatori di accuratezza

Esistono diversi indicatori di accuratezza predittiva. In questa analisi verrà confrontato l’errore quadratico medio (RMSE) di entrambi i modelli. RMSE (Root Mean Squared Error, errore quadratico medio): deviazione standard dei residui che indica la loro distribuzione; un valore RMSE minore indica un risultato migliore.

Accuracy Metrics of Naive method

accuracy(fc_na) 
##                     ME     RMSE      MAE        MPE     MAPE MASE        ACF1
## Training set 0.1052288 2.084777 1.387875 0.09408151 1.469982    1 -0.06389016

** MPE

Accuracy Metrics of ARIMA

accuracy(fc_non)
##                        ME    RMSE      MAE         MPE     MAPE      MASE
## Training set 3.988944e-05 2.07685 1.379672 -0.05328961 1.461125 0.9940889
##                      ACF1
## Training set 0.0008176214

Forecasting con LSTM Memoria a breve termine a lungo termine LSTM può elaborare sia singoli dati che una sequenza, come ad esempio un video completo. Questa applicazione è per il riconoscimento vocale e il riconoscimento della scrittura a mano. Aiuta ad evitare problemi legati alla dipendenza a lungo termine. Il loro uso più comune è lo sviluppo del processo di apprendimento di enormi problemi. Anche la memoria a lungo e breve termine è una rete neurale ricorrente, ma è diversa dalle altre reti. Altre reti ripetono il modulo ogni volta che l’input riceve nuove informazioni. Tuttavia, LSTM ricorderà il problema per un tempo più lungo e ha una struttura a catena per ripetere il modulo. Esse interagiscono in un metodo speciale e contengono quattro strati di rete neurale. Se vuoi saperne di più trovi il link qui –> https://datascience.eu/it/apprendimento-automatico/comprensione-delle-reti-lstm/

price_apple_2018 <- getSymbols("AAPL", auto.assign=FALSE, from="2018-01-01", to="2022-12-31")
APPLE_log_returns <- price_apple_2018 %>% Ad() %>% dailyReturn(type = "log")

Dal momento che stiamo per fare previsione con un modello LSTM Per prevedere il prezzo dei giorni futuri, creiamo nuove variabili da aggiungere al nostro dataset, cosi da avere più informazioni e fare un buon allenamento

se vuoi sapere di più sul significato di queste variabili aggiuntive puoi visitare questo sito–>https://juliahub.com/ui/Packages/Indicators/0sN6c/0.7.0 oppure visitare il sito tel pacchetto usato “TTR” –> https://cran.r-project.org/web/packages/TTR/TTR.pdf

#Calcoliamo le medie mobili a 10-20-60 giorni
require(TTR)
price_apple_2018$Avg_volume_10 <- SMA(price_apple_2018$AAPL.Volume, n = 10)
price_apple_2018$Avg_volume_20 <- SMA(price_apple_2018$AAPL.Volume, n = 20)
price_apple_2018$Avg_volume_60 <- SMA(price_apple_2018$AAPL.Volume, n = 60)

# Calcoliamo la % del volume medio dei giorni sopra indicati
price_apple_2018$Volume_perc_avg_10 <- (price_apple_2018$AAPL.Volume/price_apple_2018$Avg_volume_10) * 100
price_apple_2018$Volume_perc_avg_20 <- (price_apple_2018$AAPL.Volume/price_apple_2018$Avg_volume_20) * 100
price_apple_2018$Volume_perc_avg_60 <- (price_apple_2018$AAPL.Volume/price_apple_2018$Avg_volume_60) * 100

# Calcoliamo l'intervallo tra massimo e minimo
price_apple_2018$Range <- price_apple_2018$AAPL.High - price_apple_2018$AAPL.Low



#ATTENZIONE ---- se avete problemi con questa funzione potrebbe esere per via del pacchetto dplyr... provate ad usare il pacchetto stats di R ( QUINDI USARE stats::lag INVECE DI log)

price_apple_2018$perc_change_closing <- (price_apple_2018$AAPL.Close - stats::lag(price_apple_2018$AAPL.Close))/stats::lag(price_apple_2018$AAPL.Close) * 100 # 

# Intervallo tra il prezzo di chiusura dei giorni precedenti e il prezzo di chiusura di oggi
price_apple_2018$change_from_yest <- price_apple_2018$AAPL.Close - stats::lag(price_apple_2018$AAPL.Close) # anche qui stesso problema di conflitto con dplyr (aggiungere stats:: )

# Ora calcoliamo di nuovo le varie medie mobili (MA) per il range 
# negli ultimi 10, 20 , 60 giorni
price_apple_2018$moving_avg_10 <- SMA(price_apple_2018$Range, n = 10)
price_apple_2018$moving_avg_20 <- SMA(price_apple_2018$Range, n = 20)
price_apple_2018$moving_avg_60 <- SMA(price_apple_2018$Range, n = 60)

# Calcoliamo la % dell'intervallo medio dei giorni sopra indicati
price_apple_2018$perc_moving_avg_10 <- (price_apple_2018$Range/price_apple_2018$moving_avg_10) * 100
price_apple_2018$perc_moving_avg_20 <- (price_apple_2018$Range/price_apple_2018$moving_avg_20) * 100
price_apple_2018$perc_moving_avg_60 <- (price_apple_2018$Range/price_apple_2018$moving_avg_60) * 100

# L'importo totale di denaro scambiato moltiplicato per il volume (in dollari)
price_apple_2018$cash_tradet <- price_apple_2018$AAPL.Close * price_apple_2018$AAPL.Volume

# Il volume medio di contante scambiato per gli stessi periodi di cui sopra
price_apple_2018$avg_cash_trated_10 <- SMA(price_apple_2018$cash_tradet, n = 10)
price_apple_2018$avg_cash_trated_20 <- SMA(price_apple_2018$cash_tradet, n = 20)
price_apple_2018$avg_cash_trated_60 <- SMA(price_apple_2018$cash_tradet, n = 60)

# La % del volume medio ad oggi.
price_apple_2018$Avg_Dollar_volume_pct_10 <- (price_apple_2018$cash_tradet/price_apple_2018$avg_cash_trated_10) * 100
price_apple_2018$Avg_Dollar_volume_pct_20 <- (price_apple_2018$cash_tradet/price_apple_2018$avg_cash_trated_20) * 100
price_apple_2018$Avg_Dollar_volume_pct_60 <- (price_apple_2018$cash_tradet/price_apple_2018$avg_cash_trated_60) * 100

# apertura - chiusura
price_apple_2018$nightgap <- price_apple_2018$AAPL.Open - stats::lag(price_apple_2018$AAPL.Close)

# Il Gap % di vincita o perdita rispetto ai prezzi di chiusura di ieri
price_apple_2018$night_gap_perc <- (price_apple_2018$AAPL.Open - stats::lag(price_apple_2018$AAPL.Close))/ stats::lag(price_apple_2018$AAPL.Close) * 100
price_apple_2018$perc_range_previous = abs((price_apple_2018$AAPL.Close - price_apple_2018$AAPL.Open)/(price_apple_2018$AAPL.High - price_apple_2018$AAPL.Low) * 100)
price_apple_2018$perc_range_atpr = (price_apple_2018$Range/price_apple_2018$AAPL.Close) * 100
price_apple_2018$perc_range_williams = (price_apple_2018$AAPL.High - price_apple_2018$AAPL.Close)/(price_apple_2018$AAPL.High - price_apple_2018$AAPL.Low) * 100
# Intervallo di calcolo per 1 mese

un_mese_range_perc <- rollapply(price_apple_2018$AAPL.High, 20, max) - rollapply(price_apple_2018$AAPL.Low, 20, max)
price_apple_2018$un_mese_range_perc = (price_apple_2018$AAPL.Close - price_apple_2018$AAPL.Low)/un_mese_range_perc * 100

#Le medie mobili uniformano i dati sui prezzi per formare un indicatore di trend following.

price_apple_2018$EMA10 <- EMA(price_apple_2018$AAPL.Low, n = 10)
price_apple_2018$EMA20 <- EMA(price_apple_2018$AAPL.Low, n = 20)
# Media mobile ponderata
price_apple_2018$EMA60 <- EMA(price_apple_2018$AAPL.Low, n = 60)
#La doppia media mobile esponenziale è una misura dell'andamento di un titolo
# media
price_apple_2018$WMA10 <- WMA(price_apple_2018$AAPL.Low, n = 10)
# L'EVWMA utilizza il volume per dichiarare il periodo dell'MA.
price_apple_2018$EVWMA10 <- EVWMA(price_apple_2018$AAPL.Low, price_apple_2018$AAPL.Volume)
# Zero Lag Exponential Moving Average (ZLEMA) Come nel caso del doppio
price_apple_2018$ZLEMA10 <- ZLEMA(price_apple_2018$AAPL.Low, n = 10)
# Prezzo medio ponderato per il volume (VWAP-"Volume-weighted average price") e media ponderata per il volume mobile
# prezzo
price_apple_2018$VWAP10 <- VWAP(price_apple_2018$AAPL.Low, price_apple_2018$AAPL.Volume)
# La media mobile di Hull (HMA), sviluppata da Alan Hull, è estremamente
# veloce
price_apple_2018$HMA10 <- HMA(price_apple_2018$AAPL.Low, n = 20)
# La media mobile ALMA utilizza la curva della distribuzione Normale (Gauss).
# quale
price_apple_2018$ALMA10 <- ALMA(price_apple_2018$AAPL.Low, n = 9, offset = 0.85, sigma = 6)

# salviamo il file in csv 
write.csv(price_apple_2018, file = "price_apple_2018_Federico_Guzzo.csv", row.names = F)

Test di Dickey-Fuller

adf.test(price_apple_2018$AAPL.Adjusted)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  price_apple_2018$AAPL.Adjusted
## Dickey-Fuller = -1.5744, Lag order = 10, p-value = 0.7585
## alternative hypothesis: stationary

Un modo per testare se una serie temporale è stazionaria o meno, è eseguire un test Dickey-Fuller, una serie temporale stazionaria è quella la cui media, varianza e autocorrelazione sono tutte costanti nel tempo. Al contrario, una serie temporale non è stazionaria quando queste tre statistiche cambiano nel tempo. Nel nostro caso abbiamo come risultato il p-value=0.7585, quindi NON possiamo, rifiutare l’ipotesi nulla perché il valore p non è inferiore a 0,05, questo indica che la nostra serie temporale non è stazionaria.

# Rimuoviamo le variabili altamente correlate per evitare l'overfitting dei modelli
del <- cor(AAPL_lm)
del[upper.tri(del)] <- 0
diag(del) <- 0
AAPL_lm <- AAPL_lm[, !apply(del, 2, function(x) any(x > 0.9))]

** Creiamo Train e test dal nostro dataset, vogliamo provare a fare una previsione temporale di 7 giorni tenendo in considerazione i valori predetti da quelli reali**

giorni_prev = 7
n = giorni_prev + 1
X_train = AAPL_lm[1:(nrow(AAPL_lm) - (n - 1)), -17]
# la nostra variabile dipendete è AAPL.Adjusted
y_train = AAPL_lm[n:nrow(AAPL_lm), 17]
X_test = AAPL_lm[((nrow(AAPL_lm) - (n - 2)):nrow(AAPL_lm)), -17]


require(quantmod)
# Creiamo il test di validazione dei prezzi reali dei prossimi 7 giorni 
AAPL2 = getSymbols("AAPL", from = "2022-12-21", to = "2023-01-01", auto.assign = FALSE)
nostre_date <- time(AAPL2)
y_test <- as.numeric(AAPL2$AAPL.Adjusted)

train <- cbind(X_train, y_train)
# check the number of features
dim(X_train)
## [1] 1193   16
dim(X_test)
## [1]  7 16
nostre_date
## [1] "2022-12-21" "2022-12-22" "2022-12-23" "2022-12-27" "2022-12-28"
## [6] "2022-12-29" "2022-12-30"

#KERAS Deep Learning: Backend TensorFlow Applica una rete di deep learning da strati strettamente accoppiati di uno stack lineare.

gc()
##           used  (Mb) gc trigger  (Mb) limit (Mb) max used  (Mb)
## Ncells 3583271 191.4    6235594 333.1         NA  6235594 333.1
## Vcells 6815723  52.0   12259901  93.6      16384 12259901  93.6
library(tensorflow)
## 
## Attaching package: 'tensorflow'
## The following object is masked from 'package:caret':
## 
##     train
library(keras)
## 
## Attaching package: 'keras'
## The following object is masked from 'package:depmixS4':
## 
##     fit
# ##### NOTA BENE (SE HAI UN MAC CON CPU SILICON GLI ULTIMI AGIONAMENTI DI TENSOR E KERAS NON FUNZIONANO)--> Ho dovuto eseguire il downgrade di tensorflow-macos alla versione 2.9.0 e tensorflow-metal alla 0.5.0 per eliminare tali errori. Altri su questo sito hanno segnalato problemi simili con versioni più recenti. Apple deve ancora fornire versioni funzionanti per le nuove versioni di TF.
# #https://developer.apple.com/forums/thread/722361
ker = ncol(X_train)

keras_model <- keras::keras_model_sequential()
keras_model %>% 
  layer_dense(units = 60, activation = 'relu', input_shape = ker) %>% 
  layer_dropout(rate = 0.2) %>% #We apply dropout  to prevent overfitting
  layer_dense(units = 50, activation = 'relu') %>%
  layer_dropout(rate = 0.2) %>%
  layer_dense(units = 1, activation = 'linear')

keras_model %>% compile(optimizer = "rmsprop", loss = "mse", metrics = "mse")

keras_history <- keras_model %>% fit(X_train, y_train, epochs = 100, batch_size = 28, 
                                     validation_split = 0.1, callbacks = callback_tensorboard("logs/run_a"))

link al video dell’allenamento del nostro modello https://www.youtube.com/watch?v=pmUZUtPUQYc

keras_pred <- keras_model %>% predict(X_test, batch_size = 28)
real_VS_pred <- data.frame(keras_pred, y_test)
ok <- as.data.frame(nostre_date)
df <- cbind(ok,real_VS_pred)
plot(keras_history)

gc()
##           used  (Mb) gc trigger  (Mb) limit (Mb) max used  (Mb)
## Ncells 3780082 201.9    6235594 333.1         NA  6235594 333.1
## Vcells 7185229  54.9   14791881 112.9      16384 12259901  93.6
real_VS_pred <- data.frame(keras_pred, y_test)

colnames(df) <- c("DATE","KERAS PRED", "REAL PRICES")
# 
p4 <- ggplot(real_VS_pred, aes(nostre_date)) + 
  geom_line(aes(y = keras_pred, colour = "keras_pred")) +
  geom_line(aes(y = y_test, colour = "real_prices")) +
  geom_point(aes(y = keras_pred,colour = "keras_pred"),size = 2) + 
  geom_point(aes(y = y_test, colour = "real_prices"),size = 2) +
  labs(title = "Keras (Predicted vs Actual)", x = "Date", y = "Daimler Share Price in $") +
  theme_solarized(light = FALSE)+
  dark_theme_gray() 
## Inverted geom defaults of fill and color/colour.
## To change them back, use invert_geom_defaults().
p4

gc()
##           used  (Mb) gc trigger  (Mb) limit (Mb) max used  (Mb)
## Ncells 3842144 205.2    6235594 333.1         NA  6235594 333.1
## Vcells 7286576  55.6   14791881 112.9      16384 12259901  93.6
require(kableExtra)
kable(df) %>% kable_material_dark(bootstrap_options = "bordered", full_width = F, 
                                      position = "center") %>% column_spec(1, bold = T, color = "red")
DATE KERAS PRED REAL PRICES
2022-12-21 135.2805 135.45
2022-12-22 132.0124 132.23
2022-12-23 137.5526 131.86
2022-12-27 133.9008 130.03
2022-12-28 131.3170 126.04
2022-12-29 136.2456 129.61
2022-12-30 136.1113 129.93

CONCLUSIONI

I dati delle serie temporali sono ovunque e la maggior parte di essi può influenzare la nostra decisione. La scomposizione dei dati delle serie temporali può darci una visione più dettagliata del modello dei nostri dati. Analizzarlo ci aiuterà a prendere decisioni sui nostri dati. Tuttavia, la decomposizione in R può essere fatta solo all’oggetto ts dove la frequenza verrà specificata. Dopo aver analizzato le nostre time-series, di solito vogliamo prevederli per aiutarci a prendere decisioni sul futuro. Tuttavia, bisogna tener sempre presente che il risultato delle previsioni non può mai essere accurato al 100%. La previsione può essere fatta usando molti modelli e algoritmi. In questa analisi, abbiamo usato il metodo NAIVE, il modello ARIMA, HMM, LSTM. Nella nostra analisi Il modello migliore risulta esser LSTM, in quanto è quello che si avvicina di più ai dati reali, sfortunatamente, il nostro modello non è ancora sufficiente per prevedere con successo le serie temporali del mercato.