Forecast with the Adanic Method

 Abstract

In the present study, we expose a new technique of quantitative forecast called Adanique. First of its kind, because it wants to be a newcomer in the classic prediction methods known until today. Generally divided into two, in this case statistical methods (endogenous approach) and econometric methods using a purely exogenous approach. The new so-called adanic forecasting technique uses this double principle to allow a new forecasting approach. The approach of the adanic method and its demonstration on the scientific values are spread out in the sections below.

Keywords: Forecasting, Econometrics, Statistics and Adanic

Conclusion

The adanic forecasting method is both based on the ethical principle of statistical and econometric forecasting. The future estimation of the values of a chronicle indirectly by explanatory variables is the main part of the econometric forecast. With regard to statistics, we are satisfied only with the historical values of the series in question. As a result, indirect and direct prediction is the principle of adanic forecasting. Experience has shown that the result is positive because it has shown a retractability of historical and observed values of the series over the years. On the other hand, it should be recalled that the forecasting indicator did not prove very favorable to adanic forecasting. It must not be thought that there is a magnitude whose reliability is absolute, because the future behavior of a series is primarily random and not deterministic. Can it be that the average of three methods or two methods is also an option for predicting a given magnitude?

https://www.researchgate.net

Abdi-Basid ADAN

The Abdi-Basid Courses Institute

Modèle à Correction d’Erreur à Equations Simultanées : en Cas Standard et Spatial des Données de Panel

 Résumé : 

À la suite de la déficience d’un modèle performant de premier plan en économétrie, le monde scientifique semble à jamais opposé sur plusieurs théories économies, parmi lesquelles la courbe de Kuznets Environnementale, l’existence des convergences conditionnelles des économie…etc. Cette crise dont souffre le monde scientifique particulièrement dans le domaine de l’économie m’a poussé à réfléchir sur un nouveau modèle économétrique plus complexe contournant plusieurs défaillances des autres modèles : c’est le modèle à équations simultanées à correction d’erreur standard et spatial des données de panel. 

Mots Clés : modèle à équations simultanées à correction d’erreur standard et spatial des données de panel, courbe de Kuznets environnementale.

Conclusion 

La réduction en une équation réduite d’un modèle à équations simultanées exigerait une cointégration afin d’en déduire un modèle à correction d’erreur à équations simultanées. Cette approche va certainement proposer des estimateurs encore mieux efficients pour permettre de trancher sur plusieurs articles contradictoires. De plus l’introduction de la matrice des poids dans un tel modèle ne fera que renforcer ses capacités dans la prise en compte des effets spatiaux. Ici, la correction d’erreur est venue avant l’intégration de la matrice des poids. Bien qu’il soit en cours de développement, une équation spatiale à correction d’erreur ou une équation à correction d’erreur spatiale ne serait pas identique. En outre, la prise en compte des aspects plus complexes dans un modèle performant palliant les défaillances des autres modèles permettra de proposer un consensus dans un monde scientifique perplexe. Le modèle par excellence serait celui à correction d’erreur à équations simultanées spatiale ou standard en panel.

https://www.researchgate.net/publication/324978703_Modele_a_Correction_d'Erreur_a_Equations_Simultanees_Cas_Standard_et_Spatial_en_Donnees_de_Panel

Qualitative Data Analysis (QDA)

 From the econometric point of view: 

On the econometric level, the analysis of the qualitative variable (categorical or nominal) takes place in two ways: either it consists of considering it and apprehending it as being an endogenous variable, or of course retaining it as an exogenous variable and to study it in a model of qualitative econometrics. Added to this is the possibility of analyzing the link between two qualitative variables through Chi-square test of independence (especially not to be confused with the Chi-square test of adequacy). In addition, models of qualitative econometrics are diverse, and some may be very complex than they are. Let us quote a few of the simplest to the most complex: the Probit binominal model, Logit binominal, Gombit, simple Tobit, Tobit generalized type I, II, III, IV, V, Tobit double censored, Tobit multiple censorship (truncated or limited), Heckit ... etc. 

The data structure, although dynamic on the temporal component, may vary according to each observation and in this case takes the form of panel data. We will explain in the previous section some principles of these models. To find out about the nature of the model that one has to do, one can easily define according to the domain of definition of the function but also of the modalities provided by the endogenous qualitative variable. These are variables, which are often derived from inquiries concerning a character of appreciation, opinion, satisfaction. ... etc. Others can be calculated and not observed as such. Their studies are as indispensable as the quantitative ones to make significant efforts in a specific direction. 

The simultaneity of a qualitative endogenous variable explained by another qualitative exogenous variable is certainly possible. To study in such a case, it will be necessary to consider one of the modalities of the variable while retaining the others to serve as references in the meaning of the results obtained. As for the qualitative endogenous variable, it is presented in terms of probability. To achieve this, a latent continuous variable is enough to facilitate the calculation in terms of probability of the modalities of the endogenous variable. A positive probability is indeed synonymous with a growing number of chances. It is eminent to note, on the other hand, the convergence of the solution after iterations in both the concave and convex cases. Beyond the significance by predictor variable, the overall significance or the adequacy of the model is a notion of appreciation of the model design. Thus, the coefficient of determination indicator of Mc Fadden and Hekman also called pseudo Rsquare makes a judgment on the quality of the fit of the model. In other words, the explanatory power or the part of the fluctuation explained by the variables retained in the model. 

On the other hand, the remaining percentage is usually less than 50% and corresponds to the relevant variables not considered. The Hosmer Lemeshow test also goes in the same direction on the quality of the fit. Indeed, it is by the marginal effects that we know more about the impact of each variable introduced into the model. The obtained estimators give an idea of the nature of the influence of the exogenous variable on the endogenous. In principle, it is recognized that there is a difficulty in interpreting the terms of the explanatory variable. It is then better to set a modality as a reference and to interpret by comparing with those used as references in the analysis. In most cases, the choice of modeling between Logit, Gombit (or extreme value) and Probit is done with the predictive power of the model. The best of them will be retained for the final modeling. 

  From a statistical point of view: 

In statistics, there are a variety of univariate and multivariate analysis procedures, including the family of factorial methods: factorial correspondence analysis, multiple correspondence analysis, multiple factor analysis, factorial analysis of mixed data. In some of these methods, it could concern both the qualitative and quantitative variables, so we evoke a mixed analysis of the variables. In addition, other techniques study the link between nominal and mixed variables. In this case, for example, the Chi-square test, the Cramer Rao coefficient, the correlation ratio, the analysis of the variance (ANOVA), etc. Statistically analyzing the qualitative variables is to perform the same operation seen in econometrics. In other words, make each modality a new variable. Therefore, in some studies we talk about the existence of a complete disjunctive table or Burt's chart. In addition, the match name refers to the link between the nominal variables. The search for axes expressing more meaning to the data is the common denominator of all factorial methods. It is rather in Multiple Correspondence Analysis that there is a massive loss of information and therefore the need to take certain results with care. The principles of the factorial analysis of correspondences, as the name indicates it allows to highlight the correspondences between two qualitative variables. In other words, the link whose modalities intervene and specially to identify the nature of the link that can be attractive, repulsive or independent. In this sense, it is an exploratory method, descriptive of data, established by Benzecri in the 70's. 

The idea is to translate the proximity of modalities as a link between the variables and specially to grasp as an identical profile for individuals which they describe. In contrast, the multiple analysis of the correspondences is a generalization of the factorial analysis of the correspondences, which itself is a double Analysis in Principal Components on on the one hand the line profile and on the other hand the profile column in a table contingency. Another aspect of distinction is that in simple factor analysis, the raw table is not studied directly, that might be interpreted as differences between rows and columns. It is also important, when interpreting, to avoid marginal low-margin modalities lest it influence the contributions of others. In discriminant analysis, it requires the presence of a qualitative variable with several quantitative variables. 

The principle is to put in place a linear combination of the quantitative variables separating at best the studied population. The discriminant function can be obtained using multiple linear regression. According to a threshold and the modalities of the qualitative variable, one determines the points individuals misplaced. The ideal is that it requires more than a minimum number of misplaced individuals. For this, one should think about including other variables in the regression and repeating at several iterative. The particularity of the discriminant analysis is that it is apart from its exploratory function, a decision-making method. 

 From the point of view of hierarchical classification: 

From the point of view of the ascending or descending hierarchical classification, the implementation is made possible by multitudes of metric distance calculation algorithms, among which the Manhattan method, weighted distance, Ward ... etc. The idea is to reduce the number of classes by iteration by grouping the one that is similar or the one whose dissimilarity is minimal (according to the index of aggregation). In other words, we try to minimize intra class variance. This partitioning derives from the distance matrix in a space of R power variable number. Identical profiles show individuals with the same preference for a given choice or the same profile for a characteristic of the individuals sought.

Abdi-Basid ADAN

 

The concept of Darwin, a partial and partial theory

Charles Darwin, palaeontology and English naturalist, published in 1859, the book "From the origins of species." Although independent thinker, he had to borrow the hypothesis emitted by Jean-Baptiste de Lamarck, 50 years rather. In fact, almost all of his research is based on three major points: unity, species diversity and natural selection. In truth, the idea that Darwin puts forward his argument is simply illusory, because indeed, in his formulation, he suggested that one constantly keeps the excavations of remains of spaces to support his postulates.

It is very damaging to see the birth of a Darwinism current that persists with determination, educating the new generations of values ​​that do not conform to the principles of scientificity and that defends truthful facts as mere facts. In spite of what, I would demonstrate in this article the unfoundedness of Darwin's theories on the human lineage in general.

From the point of view of the uniqueness of human beings, it is necessary to ascertain by nature the true existence of a relative similarity between several living beings more on the anatomical than on the morphological level; chromosomal; molecular ... etc. As well for the species of the same medium as those of the distinct environments. This should not, under any circumstances let us think of the idea of ​​a common ancestor. The evolution of a being "origin" or plural as presupposes Darwinism can not explain the settlement of the billions of billions of living species on Earth. This constitutes a misinterpretation of nature and its evolution. In principle, the commonality between species is irrefutable and must be seen from the angle of belonging to the same supreme being and not to the same common ancestor. Now it is amazing to explain an infinite multiplicity of beings, knowing that others are undiscovered until 2018 as being a drift of one or more common ancestors. On the principle of science, phylogeny is not in line with the evolution of species in a much more global dimension, it is partial and partial in its concept. Time would not be the mastermind of the "magic" differentiation of these species. In the excavation of skeletal remains, discoveries constantly call into question even the relative consideration of Darwin's principles. In 2004, a team of prospecting Aramco Exploration Team would have discovered, besides in Saudi Arabia a human skeleton of outsized proportions. The finding of remains of the giant men would be the incontestable proof of the basis of Darwinism thought.

With regard to the diversity of living things, it is quite clear that nature is balanced by this principle, it is necessary and unconditional diversity to allow the longevity of life on Earth. In this sense, all species intertwine with each other. This diversity can not be denounced, it is welcomed without exception, Again, it must be remembered that plants are living beings and that the earth is so. I can not understand how Darwin could think of these beings. Descartes in his explanation brings about the origin of the mistakes that rushing further pushes people to make mistakes. Darwin's attitude is the perfect example. This diversity, in fact, derives from the ingenuity of a supreme being unequivocal and unambiguous. Finally, natural selection, although defined as the equilibrium of the environmental balance over species over time. It is quite complex to be able to understand in a complete dimension. Reproduction of some species is in favor of others. The disparity of others is advantageous for future new species, which can derive from minority strains. The complex environment, predator and prey governs this complex principle par excellence.

In conclusion, to think that uniqueness gives way to the diversity of species is illusory. Although natural selection is the driving force behind the evolution of species over time. It is not the latter that determines the period between the uniqueness and the diversity of species. The mystery of the uniqueness of living beings still remains to be seen, in my opinion it is the manifestation of a supreme being which defines it. The diversity of beings is an essential essence for the stability and continuity of life on Earth. The role of natural selection should be understood as the interaction between environment and species. From these points of view, there is no common ancestor in the human lineage with apes. The only similarity can not justify this, it is the utopian drift of the currents Darwinisms.

 The perception of the notion of time

A prodigious question does time really exist? the reality in front of us, would it hide a truth as intimate as it is trying to conceal it from the good of a revolutionary discovery of all time. Let's go back to the 17th century, from one evening in 1666, under a tree with a clear sky, Newton observes the fall of an apple, so why would not the moon fall like that? was born classical physics. It is remarkable and amazing to see that nature in the evening of 1666 transmitted its secret to Newton. Even further, in 1905, Einstein extends the classical physics to the infinitely small, at the very heart of matter with a limited relativity (in a frame of reference not subjected to a fictional acceleration as the non-Galilean case) and on a cosmic scale with its general relativity (existence of the curvature space temp). Here, if we see that we are still in the explanatory trend of nature in two different dimensions. The most disturbing question for Einstein was the light and it is from this interrogation that the door of his new physics opens, quantum physics. In truth, he was the only one to understand that time could slow down from one point to another. His explanation of the photoelectric effect or the absorption of photon (quantum of light) by the material, when interacting with the light, earned him the Nobel Prize in physics in 1921. An extraordinary physics, which is not more superficial as it was with Newton but leaving, this time in depth of nature to explain the wonderful laws that govern. From this point of view, it is clear and unequivocal that science in general is the alternative for man to better understand the nature that conditions every element of life. This is undoubtedly the key to better understanding the main condition of life. In this sense, there is no science without nature.

What does time depend on? The paradox that I think one can emit is to say from the revolution of the earth and that of rotation, the count in 24 hours and the 365 days of the year is concrete but separately. If the 24 hours and 365 days are not excesses of the same movement, we have a discordance purely clear between the short term and the long term. In this sense, reference or position and mobility play an important role in determining time. It is therefore very relative and not absolute.

 

 

To discover the secret on the principle of time, we will first have to lend an intention to nature. Because the latter can serve us a third door on the mystery of physics. Time, by definition, is nothing more than the interaction between life and nature through the environment. It requires an environment conducive to life under the order of nature. Perfectly knowing the non-continuity of life, all the stakeholders of life are recycled and do not evolve on an infinite trajectory. It is just renovating in a new dimension exactly as the bell curve of Laplace Gauss says. Some species live longer than others. Time is not only unevenly distributed over every point of reference, but there are also parts where the time is even more abnormal. Why this inequality or super inequality of time in space? One could consider the idea that nature although it is based on the principle of balance, adapt the time to maintain its foundations (stability, diversity cycle, balance ... etc.). 

Obviously, time does not exist for certain entities. It is a primary component to which time does not apply to it. This idea is indebted because the element in question is not related to nature or life. It is undoubtedly in this case the vital breath in other words the soul. To survive, he has no need to drink or eat. Indeed, it is the body that determines the need but not in him. He is therefore outside the principle of life. Nature does not apply to him. And from this point of view, time too. Many physicists ask themselves whether time does not exist. To tell the truth, time does not exist for certain entities, but it is for others in a relative way. Even more surprising, say that the grandfather of a family and his grandson have the same age spiritually. It's different, it's the age of their organisms.

In a word, time is the stranglehold of nature on life. It is relative according to the referential chosen and applies only on all the elements which contribute the life, the solar system and its various forms, the species with different life expectancies, the cycle of the environment and time and space itself are all concerned.

Spatial Autocorrelation Patterns and Inequality Analysis

Objective Description

This analysis examines spatial autocorrelation patterns in synthetic raster data using the Getis-Ord Gi* statistic to identify hotspots and coldspots, and measures spatial inequality with the Gini index. It involves generating a synthetic raster (representing NDVI), visualizing it as an NDVI map, creating a shapefile, extracting raster values, computing spatial neighbors, calculating Gi* statistics, classifying hotspots, visualizing hotspots with a thematic map, and exporting outputs as shapefiles and rasters. The analysis highlights spatial clustering and inequality in a simulated dataset.

🎯 The detailed methodology and results can be accessed through this link:

👉click here now! : https://rpubs.com/abdibasidadan/spatialautocorrelation

Abdi-Basid ADAN, 2025

The Abdi-Basid Courses Institute

Multi-Objective Optimization Analysis with NSGA-II

Objective Description

This analysis employs the Non-Dominated Sorting Genetic Algorithm II (NSGA-II), a robust multi-objective evolutionary algorithm, to address complex optimization problems. The study includes two scientifically relevant case studies:

1. Car Example: Optimizes fuel consumption and maximum speed based on vehicle weight and power. This is critical in automotive engineering for designing sustainable vehicles, balancing environmental impact (fuel efficiency) with performance (speed), which influences market competitiveness and regulatory compliance.

2. DRASTIC Index Example: Optimizes weights of the DRASTIC parameters (Depth to water, net Recharge, Aquifer media, Soil media, Topography, Impact of vadose zone, and Conductivity) to maximize correlation with nitrate concentration (NO3) while minimizing Root Mean Square Error (RMSE). This enhances groundwater vulnerability assessment, providing valuable insights for environmental management and policy-making in regions prone to contamination.

The analysis generates Pareto fronts to visualize trade-offs, computes optimal solutions, and exports results as CSV files for further scientific evaluation.

https://rpubs.com/abdibasidadan

The Abdi-Basid Courses Institute

Hierarchical Cluster Analysis with Dendrogram

Objective Description

This analysis applies hierarchical clustering to a simulated dataset of 50 observations with 5 variables using Ward’s method and Euclidean distance, visualizes the results with a dendrogram, and extracts cluster assignments.

https://rpubs.com/abdibasidadan/Dendrogram

The Abdi-Basid Courses Institute

Monte Carlo Simulation and Bootstrap for Estimating the Cost of Wind Energy

Objective Description

This analysis estimates the Levelized Cost of Wind Energy (LCOE) using Monte Carlo simulation to account for uncertainties in key parameters (investment cost, operation and maintenance cost, interest rate, capacity factor, and lifetime) and bootstrap methods to compute confidence intervals for the mean LCOE. It includes wind speed modeling, vertical extrapolation, power density calculation, Weibull distribution fitting, and correlation analysis.

https://rpubs.com/abdibasidadan

The Abdi-Basid Courses Institute

 Unscientific science

The most beautiful expression of the day would be to hear and expose itself under a denomination of absolute confusion in the foreground, however, towards a purification of truth to take precedence throughout this work of reflection. It is unequivocally undertaken with the promulgation of a neuralgic point of science, strictly speaking the unscientific science. The underlying idea that one would like to imply par excellence is much more complex than the perplexity it raises in us with simple reading. With good courage and clarity of preference, the reflection contained in this article derives from an assertion which will leave us in no way, without any indifference. For the sake of simplicity and introspection, paying attention to relationships that may well exist between three key entities that are formally daily. In this case they are "the human being", "life" and "reality". The chosen research goes in the direction of first to diverge the perspectives in order to converge them all towards a common thesis, which only deals explicitly with the heading of this present report. In any case, essential is it? that our supreme convictions, at a certain moment, at a certain stage of life, should lead us to detach ourselves from the height of the grandeur of consideration as a fact which is true and without absolute contradiction. Irrational to believe that the human being really exists in nature. Indeed, the recognition of the human species is a question of form and not of substance. If we come to distinguish an elephant from an ant, it is because, first, the principle of morphology is fundamentally based on this notion of classification. Very often, any object in a room is deprived of the specific character of beings: that of self-determination. In other words, we all agree that a gust of wind and eye contact with another person are not interpreted in the same way by the mind. There is always this agreement of priority to the thinking entity of a value which itself resituated by its faculty of self-determination, a kind of prestige by nature. But it is not through a soul that we can afford. A lifeless body is not so different from any object. To signify here, that the human being is an invisible species (soul) that seeks a visible identity with this body which serves him to prove his existence on Earth. As a result, the human species exists only by and indirectly through physical visibility. 2 The question of life is the same allegation. I do not attach much here to the subjectivity of the principles, but privilege more and unconditionally the universality of the latter. It makes sense that we are continually caught between days and nights throughout the Earth's journey. You have never wondered if this is the case. To tell the truth, the darkness of the universe makes us plunge into the night and the sunlight into the day. There is no day and night proper, but it is rather the ambush or "dance" of the Earth that gives rise to this sensation and that is so. To simply say that life is this illusion that suggests that everything is normal, logical and real. In this sense, life itself does not hold to logic from several angles. From this condition, does reality constantly mock our consciences? Not at all, she is known to be stable, constant. So where does the disturbing reality come from? From there, starts the cogitation on unscientific science. What is science first and foremost? It is the ability to interact with nature, to understand better, to grasp it, above all else, to want to domesticate it in several forms at will and to the satisfaction of man. The quest for knowledge began the day when man knew the heavy burden of life on Earth. In these precise terms, science is the one that makes us discover reality. But if the latter is utopian, in this case, did not it be wrong in advance on this basis? Quite simple, how to explain with science some plausible eventualities? Like for example a land passage in the Atlantic Ocean as it was in the Red Sea during the exodus of the Hebrew people; the steps of the mountains; the compatibility of incompatible natures; the birth of a sheep by a cave ... etc. It cannot be said that it belongs to the divine manifestation. Absolutely, although we are constantly explaining, we must note that the reality is not as real as it is. Many things are hidden from us. So much so that everything is done to evaluate us and that a much different story exists beyond reality, the human being and life. It should be remembered at the end that science, although an admirer of nature as such, it can sometimes prove to be sometimes contradictory to different degrees. As we have just demonstrated, all that we know of reality is far from being "a logical whole". Once again, certain natural phenomena already discovered persist in the most total blur, this is the case of the rock of Jerusalem is suspended in the air, the mysterious holes that form on the surface of the earth of recent times ... etc. 

Abdi-Basid ADAN

 Integration and Cointegration of Variables (ICV)

Abstract 

In this report, we present an overview of notions of stationarity (integration) and cointegration of variables. Although they occupy an important place in statistics and econometrics, their concepts have, for many researchers, been a little vague. I allow myself to further substantiate this by proposing only the essentials from the theoretical and empirical point of view on the integration and cointegration of variables. 

Keys Words: ARDL, Box-Jenkins, Engle & Granger (1987), Johansen (1988)

 Introduction

 A temporal variable, a process of short or long memory requires, before being analyzed, a study making it possible to raise the great characteristics of a quantity from the statistical point of view. This is among others, its tendency; seasonality; its stationarity, its law of probability or its density, etc. One of the most interesting and in most cases very interesting is the question of the stationarity of the chronicle. Several processes are at our disposal for its implementation. In addition, the interest of cointegrated variables in the case of a linear combination has become viral in econometric studies. For fear of presenting an illusory model with no-standard coefficients, the literature hurries to propose a lot of tools to come around this major obstacle in econometrics. In the following sections, we present the stationarity of the variable in a more concise manner and discuss the notion of cointegration later. 

Abdi-Basid ADAN * Integration and Cointegration of Variables (ICV) * abdi-basid@outlook.com P a g e 4 | 8 I. 

Integration of Variables (IV) One of the fundamental concepts in statistics and econometrics is actually that related to the integration of variables. What is really the integration of variables? It turns out that in statistics, the idea of no-stationarity is rather commonly used compared to that of integration of variables. However, the meaning of the latter is versatile from one discipline to another. In statistics, a temporal process fluctuating independently of time is said to be stationary. For illustration (chronogram and correlogram), we can witness a concentration of the chronic around its mean value (mathematical expectation) or a rapid decrease of the autocorrelation function. In terms of probability, the joint distribution (law or density) is identical for the first k and k + 1 variables. In this sense, the process is stationary in the strict sense (also called first-order or strong). To say more, whatever the moment t considered, the variation of the chronic in t and t + 1 is not influenced by the temporal reference. From this point of view, its properties, also known as "stationarity conditions", namely mean, variance and covariance, are all convergent and independent of time. In view of the complexity of estimating a law of probabilistic nature for a given distribution, it has been proposed to remedy this, a second sense of stationarity, called second-order stationarity (weak sense or covariance). Only the conditions of expectancy and the variance of the variable are required and indispensable to judge and affirm stationarity against a variable. Knowing the third condition also includes the second when k is zero (cov (xt, xt-k)). This ensures a pioneering role especially in the forecast of a time series. In principle, the confidence interval of the estimated values depends on its validity. It is mentioned in the literature the existence of a multitude of forms of non-stationarity of a temporal sequence among which the stationarity with tendency (TS: Trend Stationnarity, Nelson et Plosser (1982)), which graphically, one witnesses a growth evolution of the series during of time: as a result, stationarity conditions are not met. To circumvent this adversity and later obtain the series devoid of trend component, we have to choose one of three techniques including straight ordinary least squares, simple moving average or by the filter Hodrick-Prescott. Here we obtain a white noise Hamilton (1994), by definition 

Abdi-Basid ADAN * Integration and Cointegration of Variables (ICV) * abdi-basid@outlook.com P a g e 5 | 8

 stationary but of second order. In the statistical vocabulary, one also evokes the independent white noise, because of the independence between the variables, which implies the nullity of the value of the covariance or simply the decorrelation between the variables (the reciprocal being already false). Sometimes of a Gaussian nature, it is both independent and follows a normal law. Under these conditions, the white noise is strictly stationary (or stationary in the strict sense). In contrast, a stochastic process is not necessarily stationary. The second form of stationarity, for its part, is of stationarity by difference (DS: Difference Stationnnarity, Nelson et Plosser (1982)). It is due, in particular, to the influence of the series by its own historical values. In evidence, the differentiated in a good order differ from zero allows to stationise the series. In practice, two families of stationarity test are announced. It is, on the one hand, the test having the null hypothesis of stationarity (KPSS Test) and, on the other hand, the tests having the null hypothesis of non-stationarity (Dickey-Fuller or Phillips-Perron test). Even more, a variable can be stationary with constant and trend; neither one nor the other or with constant only, according to the significativities of these parameters. The nuance which it seems appropriate to specify is that of telling oneself how a series can be stationary with a tendency, whereas it is the tendancy component that must be deprived of it. Quite logical, in fact, it would mean the same thing and not seen as two different or contradictory propositions. In a word, variable integration consists of depriving the time series of any trend or historical disturbance in order to identify the stable series for use in various analyzes. Finally, it is with the support of this primordial and imperious step that we can afford to begin the estimation of the following models precisely models of Box-Jenkins (ARIMA, SARIMA, VAR, ARMAX, ... etc) or Engle (ARCH, GARCH, ... etc). 

Abdi-Basid ADAN * Integration and Cointegration of Variables (ICV) * abdi-basid@outlook.com P a g e 6 | 8 

II. Cointégration des Variables (CV) Cointegration is indispensable in two-dimensional analysis both in statistics and econometrics. We are aware of the imminent error that can lead us to the simultaneous consideration of two or more variables without first having used a study by variable. In a way, "better first evaluate each variable before the considered helpful". Introduced in economics by Engle and Newbold (1974), then developed by Granger and Engle (1987) and continued in 1991 by Johansen, the notion of cointegration is indeed, the integration of the linear combination between two variables taken together to which tests in the literature allows us to assert in case of presence of cointegrating vector. This is the Engle and Granger tests (1987); Johansen (1988); Johansen and Juselius (1990), ... etc. It is so common to end up with the problem of fallacious regression in the modeling of the variables, because of a linear regression on non-stationary variables, which in principle testifies to an explanatory mechanism R² and a t of Student very appreciable when in reality, there is no relationship between them. In this perspective, the distribution of the estimated parameters no longer follows a Student's law but a Wiener process (Brownian motion), which occurs when the variance of at least one of the variables diverges, due to its dependence on the temporal dimension and well identifiable with recursion procedure. Indeed, good prospects of the residual analysis accompany for the final validity of the model. When the residue is not stationary, it is assimilated to the case of presence of autocorrelation between the model residues. The order of integration of the residual model is not necessarily below that of the variables of the model. In obviousness, a residual component of the stationary model with an order different from 0 is considered as a model where the autocorrelation of the residues is imminent. In other words, autocorrelation and stationarity are indirectly involved through the last condition of covariance (in the weak sense of stationarity). The goal of cointegration is to be able to determine a stationary residue while working on two nonstationary variables in level. The idea proposed by Engle is in the sense that in the short term the variables diverge, but there exists in the long term a stability, a balance between them. A common evolution of variables. Now knowing the possibility of no-zero order cointegration, in the long 

Abdi-Basid ADAN * Integration and Cointegration of Variables (ICV) * abdi-basid@outlook.com P a g e 7 | 8 

term the desired perspective exists according to well-defined conditions upstream and downstream. Among them, the same order of integration of the variables, verifying the possibility of cointegrating the variables, in other words an order less than or equal to the common order of integration of the variables. In addition, downstream, a negative sign and significant equilibrium restoring force coefficient (or speed of adjustment) is required while checking the stationarity of residues obtained. Gran-ger's (1983) theorem highlights the cointegration relationship and the error-correction model. Short-term relationship estimation with OLS is only possible when the variables are differentiated. In other words, by integrating in the model, the delayed variables as explanatory. On the other hand the long term relation are es-timées in level by OLS. No method is perfect, in this sense the disadvantage of the Engle-Granger approach (1987) is the non-distinction of cointegration relation. In principle, it has only one cointegration relation, whereas it is of number k-1 relationship with k the number of variables. Johansen (1988) provides a solution to this problem with his multivariate maximum likelihood approach. For the validity of the VECM model, a cointegration rank of less than the number of variables and non-zero is required, which is evidenced by the maximization of the likelihood log. Otherwise, a VAR model (p) is estimated in place of a VECM. On the other hand, as in the self-gressive Vector model, the specification of the model according to the absence or presence of a constant and trend is necessary. We can determine them with their significativities once estimated. The one-step method of BANERJEE et al. or MCE at Handry makes it easier to interpret the long-term relationship. In addition, to estimate long-term relationships in the case of small samples, the two-step model could lead to estimation bias according to Banerjee, Dolado, Hendry and Smith (1986). The ARDL model "AutoRegressive Distributed Lag / ARDL" and the cointegration test at the terminals of Pesaran et al. (2001) take a new approach to overcome the case of cointegration of variables at different orders of sta-tionality. On the other hand, when the variable integration order is greater than 1, the application of the ARDL model poses a problem. 

Abdi-Basid ADAN * Integration and Cointegration of Variables (ICV) * abdi-basid@outlook.com P a g e 8 | 8

Conclusion 

Stationarity is more than ever a prerequisite for studying variables before introducing them into larger studies. It involves denying the temporal process of trend and historical disturbances with appropriate procedures to conduct statistical and econometric analyzes. In the same way, cointegration proves even more complex, because one wants to be interested in working on level variables, while avoiding that the regression is misleading. Several approaches are pressing, one more efficient than the other. In both cases, these notions can not be taken with negligence, because a scientific production based on the temporal process inevitably depends on it. 

Bibliography : 1. Atoumane Diagne (2015) : Econometric modeling of electricity consumption in Senegal from 1999 to 2015 2. Hélène Hamisultane: Error-correction model and applications 3. Jonas Kibala Kuma (2018) : ARDL Modeling, Cointegration Testing and Toda-Yamamoto Approach: Elements of Software Theory and Practice 4. Lonzo Lubu Gastonfils (2015) : Application of the Box-Jenkins Prediction Method

ECONOMIC GROWTH, CO2 EMISSIONS AND INCOME INEQUALITY: THE NORTH AMERICAN CASE

Resumé: 

La croissance économique, les émissions de dioxyde de carbone et l’inégalité de revenu sont les thématiques le plus controversées de nos jours auxquels nous développons davantage dans cet article. Ayant déjà réalisé le même sujet dans le cas de l’Afrique, le présent travail s’applique sur un continent riche celui de l’Amérique du Nord (Les Etats-Unis, Canada et Bermudes). En réalité, il s’agit d’étudier à la fois économétriquement et statistiquement les grandeurs définissant les trois notions évoquées précédemment dans un cadre dynamique. Un modèle à équation simultanées en données de panel est à l’ordre du jour, qui regroupe trois équations : il s’agit de l’équation de Solow, Kuznets environnementale et Kuznets sociale. Enfin, nous prévoyons l’évolution des grandeurs à la fois de la croissance économique et du développement durable avec le principe de Box-Jenkins pour les trois pays cités. Mots clés : Croissance économique, développement durable.

Conclusion : 

According to the study on the concept of sustainable development of rich countries on the continent of North America. We have, thanks to a model with simultaneous equations, proved that the control of the deterioration of the environment or income inequality over the short, medium or long terms remains illusory. The functional forms mentioned in the literature are not in line with those we have developed using the least squares estimators. On the other hand, US per capita growth will stabilize around $ 4,000 over the next 10 years, unlike Canada and Bermuda. These will continue their rhythms of a still dynamic growth (more than 7000 $ and 5000 $ respectively). CO2 emissions and income inequality on the American continent will still occupy positions that are far from reassuring, except for the United States (a slight decrease for emissions). 

https://www.researchgate.net/publication/326468307_ECONOMIC_GROWTH_CO2_EMISSIONS_AND_INCOME_INEQUALITY_THE_NORTH_AMERICAN_CASE

Unpublished Perspectives on Consciousness

Consciousness is commonly defined as the property that living beings boast about in order to assert their existence, to be animate, to be able to participate, to act and to play a role, which in their eyes seems essential in life. One characteristic among many others, standing in the foreground and which consists in lamise in immediate knowledge of the state of circumstance. In classical terms, it is acting and at the same time knowing that one is acting. As opposed to unconsciousness, which manifests itself in situations of unthinking behavior or automatic achievements. Freud, says it clearly, "the ego is not the master in his own house". Except the latter, "the id," takes relay also and identifies itself as the drive reservoir and the superego, as the instance of prohibitions and ideals. 

The subject is no longer a united entity, but a self-conflicting reality. In this article, whose primary vocation is a contribution to the notion of consciousness and its derivatives (preconsciousness, unconsciousness) known until today. With this new horizon, which promises more than ever to widen the percepts of the conscience with prospects ready to be exploited, it is among others the hyperconsciousness and the hypoconsciousness. Like Christoph Colombe in the quest for America, consciousness is certainly interpreted, indirectly when the soul is conquering the psyche system. A jouissance on multiple capacities at its disposal. But, it seems beforehand and on the same wavelength of the imperfection observed in the present life, brings without sparing consciousness to define its defects and limits, thus constituting the fundamental object of unconsciousness. Like these articulations, it emerges fiercely that consciousness becomes the coordination between the psychic system and the vital breath (soul), thus keeping it in continual action through the electrified charges of the neurons. 

It is imperative to put out of the shadows and away from confusion an emblematic subject, that concerning the species of the same intangible family. In these terms, the soul is not assimilated to a spirit. In addition, the latter branched out into several subspecies, some of which seems more hostile to humans. By their wickedness, a demonic spirit can afford to replace the soul in the colonization of the psyche to different degrees. It is not bewildering to see a form of fidelity of hearing between the soul and the spirit in this life. However, it turns out that the soul does not act alone in its decision-making, in its actions. Hyperconsciousness consists, then, of the soul's own action when it is not biased by any form of direct or indirect spiritual external influence. In contrast, to hypoconsciousness, the soul is forced to yield for one or more spirits, the management of the psyche otherwise the consciousness is placed at the service of another spiritual actor. There is reason to understand that the body, as mystical as it is, becomes a mere visual identity for these invisible beings. Especially since moments of hyperconsciousness are almost similar to the state of consciousness. 

The difference with the superego is self-judgment in a neutral way about one's own acts performed. A form of regret felt with regard to an action performed for a so-called "state of disobedience" motive of search. Contrary to what Sartre says in Being and Nothingness, "Other is, by principle, the one who looks at me and through whose eyes I become aware of what I am. " Here, by contrast, it is by no means another person who freezes his own freedom of action, but his own self reveals to what degree of prosaicism he has allowed himself to exercise. Another by amplification of the notion of consciousness, it is the order of the day of clarification for a specification between human consciousness, besides the animal, the spiritual one.

In this term, the soul, once in the body, its consciousness is immediate and is not provided with memory. That is to say, she does not need experience. Indeed, a baby with its external environment shows above all that it is dominant. How is it done without an experiment being transmitted to him, manages to control? It is by the spiritual consciousness that this questioning is explained more. In the absence of the death of a person the human consciousness stops, but the spiritual one continues. Let us go even further, in the same way as Einstein, who claimed that science without religion is lame and that religion without science is blind. In this perspective, Islam confirms the lesser knowledge that is shared with us on the soul and that it is among the angels the being that is most valued. We realize that man is seen to be superior in many contexts. Example of event: the case of anti-Semitism, racism, ethnicity, including boasting.

Abdi-Basid ADAN

Phematic

In the reasoning opposite, it is a special ease of an explanatory monologue based on a novel idea of a newcomer in the history of science. But above all, taking a stance for breathtaking attention with the following expressions including the hematic, rhematic, phreatic or thematic. Notwithstanding, the phematic the main object praise in this interview, brings us back to the question of knowing, is that a phematic, or a phematician? 

Starting from this preamble, which holds as an introductory dimension of a luminous idea, according to which a possible contribution to the future echelons seems decisive. Dear readers, as you all know, science, in the general sense, is a way for us, on which we claim to have a mastery of it is that, nature, natural balance and therefore natural laws. Is it not a tremendous job for man to develop himself into an activity of reflection leading to a career of life? But to the detriment of the multiple and complex needs that constantly lurk at us, we are more than ever forced to go beyond the painful and intrusive work of research. It breaks down the major obstacles of life, like a show in which constraints take different forms for multiple reasons. Time is changing and weariness is making us return to a comfortable lifestyle vision. That said, science since the Paleolithic; the Mesopotamians, the Babylonians, the Egyptians, the man, by his capacities to evolve, is supposed not to become a victim of the nature and not for the nature, a submission to the man. Science is at the intermediary service between man and life by slamming nature. Centuries in the centuries, innovation, discoveries and the progression of knowledge levels are climbing upward. 

The logic of any kind of knowledge, although a concept known indirectly, does, in fact, only resume the structure that consists of nature. In all simplicity, any imagination so naive or not part of an idea, a questioning of reality, which surrounds us. It was the case of the moon, why would not it have collided with the Earth? Such has become the cradle of classical and quantum physics. Compared to nature, a socalled atom represents the element of first constitution with which matter is formed. This unitary or elementary profile is proliferated, in two, then in three, until generating a block of idea in an extremely close bond in order to circumvent the confusion and the refutation. In simple terms, the origin of a cause, which itself becomes the source of another and so on, in the same way as a family tree. Noticing, one thing, through which a relationship with nature is an imminent reality about the structure of a molecule with these innumerable atomic compositions. We started from an idea to formulate a theory, a set of notions contributing to a common vision In turn, this theory is associated with others, which by principle of iteration, form a mega structure of knowledge, which is called a discipline. Incontestable, to see that we are in the horizon of the infinitely small towards the infinitely large. 

From this perception, science becomes a perfectly reflective mirror of nature. In more detail, the Earth as small compared to the sun, the sun as huge as it is in the solar system becomes, meanwhile, a tiny and invisible celestial object in front of Arcturus, who in turn becomes it before Betelgeuse ... etc. The latter still certifies the irreproachable fidelity of science to the composition of nature, its order and its equilibrium. No reason to want to remind that the physical equations have good reasons to marry nature. Can we succeed in reconstituting science, knowledge, and knowledge in a disciplinary unity, thus uniting all fields of science? This has been the work perpetrated by authors for years, which finally for today, comes to see the day. I prove the possibility of designing a unique science from all the different branches of science. 

The phematic then becomes the discipline par excellence of reference, whose expression stems from the conjunction between philosophy and mathematics synonymous with the expression mathilosophy. For phematic expert, two skills are in this sense, required in particular in philosophy and mathematics. No discipline seems to escape from these complementary doctrines to each other. Philosophy, science of formulation and cogitation and mathematical science of confirmation and demonstration. By definition of things, we will always hear, by phematic, the science of demonstrative cogitation. In a word, it means by mother of science philosophy and father of science mathematics. 

Abdi Basid ADAN

How to Diagnose Change in Precipitation and Temperature with the R program ?

Introduction

The primary purpose of this study is to demonstrate the scientific feasibility of a climate change project. To do this, it is imminent to call on a climatic database of local or satellite stations. We are going to show how to use these programs of RClimDex, RClimTool, RHtest and ClimPact under R. It is noteworthy to recall that knowledge of latitude, longitude, types of variables and their faults (missing, extreme, erroneous) are necessary for to be able to ensure the required analyzes on climate change.

I.                 RCLIMDEX

The import is done in R with an import program allowing to display a window of RCLIMDEX. With the outputs of the RCLIMDEX program under R, four documents including the climatic indices; the qualities of controls and graphics and trends will be provided to us. The indices are 27 indicators divided into precipitation, maximum temperature and minimum temperature, among which:

 Temperature indices including Min and Max temperature (SU25, TXx, TNx, TMaxMean, TNn, TMinMean, DTR,,… etc.).

 Precipitation indices including humidity and rain (PRCPTOT SDII, CDD, CWD, RX1day, RX5day, R10, R20, R95p, R99p,… etc.).

The log document provides graphs of temperature and precipitation variations, as well as other graphic illustrations in Excel format for tables and pdf for images.

The plot file groups together the figures of the climate change indices associated with their respective slopes and subjected to the Mann Kendall trend test in order to evaluate the significant trends (sign of the slope; Student significance; regression line; correlation; coefficient of determination). The last Trend document sufficiently summarizes the information in graphic format from the previous file.

I.                 RClimTool

L’importation The import is done under R with an import program allowing to display a window of RClimtool. With the outputs of the RClimTool program, three files, including homogeneity; Missing_data and Quality_Control will have to be extracted to use. The first file that emerges is the Quality_Control, for each of these variables, we study the atypical, outlier or outliers values ​​presented by our daily, monthly and annual observations. The second is that of missing value, the processing will be done separately with values ​​and without missing value. The last homogeneity document presents the various tests of normality, Mann Kendall significance, correlation between sites, with the associated graphs. It is possible to deal with missing values, because this prevents the calculation loop from functioning normally. Various deterministic methods such as statistics such as the average; the median or even probabilistic techniques (Amelia, MICE, PPCA, etc.) are used to supplement missing values ​​if they do not exceed certain thresholds of significance.

 

 

II.             RCLIMTREND

 

Plusieurs Méthodes statistiques pour les sciences du climat comme le Tests d'homogénéité absolue SNHT absolu ;  SD différent, breaks, Buishand, Pettitt, rapport de von Neumann et ratio-rank, Worsley et Craddock, tests d'homogénéité relative SNHT absolu ; break SD ; la correlation ; la recherche des valeurs aberrantes sont les diverses atout que l’on peut bénéficier avec ce programme.

III.         ClimPact

Le programme de ClimPact s’exécute de la même façon que RClimDex, les issues sont similaires cette fois-ci, une étude plus approfondie que le précèdent avec une évaluation de variation extrême. D’autant plus que nous avons comme suppléments sur le treshold (seuil). Les différents de statistiques descriptives (boxplot, histogramme.etc.). En comparant les valeurs des tests de Mann Kendall respectifs sur les indices de calculs retenus, on constate que c’est presque similaire.

IV.          RHtest

With this program, it is possible to detect point changes of the studied series and to inform us about what needs adjustment. Geoclim ; Anclim ; EdGCM ; MASH ; ProClimDB…etc

Spatial analysis tools are designed for analysis in climatology. Provides scientists and non-scientists with more decision support and facilitates filling observation gaps (precipitation and temp) using a mix of field observations, readily available satellite data and networks.

Quality Control

The preliminary study on data quality detects the extreme values ​​(tools) by realizing the statistics of maximum and minimum. The existence of missing data disrupts the analysis of the observation series for calculating indices under the RclimDex program (Zhang and Yang, 2004). A month containing more than 3 missing days is not considered in the year. An incomplete year is also private in index calculations. Statistical indicators such as the median make it possible to complete the missing values ​​when the proportion of these values ​​is negligible.

         Homogenization

Several tests make it possible to detect the homogenization of the series observed. The bivariate case like the Craddock test (1979) or the univariate case like the PETTITT test (1979). The programs of RclimDex, Rclimtool, Rclimpact (Zhang and Yang, 2004). allow the analysis of the homogenization of temperature and Precipitation series.

1.   Climate Change indices

The variability and trend of precipitation and temperature is determined with the calculation of indices under the application of RclimDex. There are 11 for the Precipitation and 16 for the temperature. The definition of each of these indices and their units are presented in the following table.

 

Table1. List of 11 ETCCDMI core Climate Indices for precipitation

Indicator	Name	Définitions
RX1day (mm)	Max 1-day precipitation amount	Monthly maximum 1-day precipitation
Rx5day (mm)	Max 5-day precipitation amount	Monthly maximum consecutive 5-day precipitation
SDII (mm/day)	Simple daily intensity index	Annual total precipitation divided by the number of wet days (defined as PRCP>=1.0mm) in the year
R10 (day)	Number of heavy precipitation days	Annual count of days when PRCP>=10mm
R20 (day)	Number of very heavy precipitation days	Annual count of days when PRCP>=20mm
Rnn (day)	Number of days above nn mm	Annual count of days when PRCP>=nn mm, nn is user defined threshold
CDD (day)	Consecutive dry days	Maximum number of consecutive days with RR<1mm
CWD (day)	Consecutive wet days	Maximum number of consecutive days with RR>=1mm
R95p (mm)	Very wet days	Annual total PRCP when RR>95th percentile
R99p (mm)	Extremely wet days	Annual total PRCP when RR>99th percentile
PRCPTOT (mm)	Annual total wet-day precipitation	Annual total PRCP in wet days (RR>=1mm)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

R Program Analysis Climat Change

 

#TREND ANANLYSIS

data=read.csv(file("clipboard"),header=T, sep="\t", row.names=1)

data

str(data)

 

#install.packages("trend")

require(trend)

 

#data(maxau)

prectot <- data[,"prectot"]

mk.test(prectot)

 

dat=read.csv(file("clipboard"),header=T, sep="\t", row.names=1)

data(nottem)

str(dat)

dat=round(dat,1)

dat1 <- ts(dat, frequency=12, start=c(1961,1))

plot.ts(dat1)

 

logsouvenirtimeseries <- log(dat1)

plot.ts(logsouvenirtimeseries)

 

require(trend)

smk.test(dat1)

 

#correlated witn precedent month

csmk.test(dat)

 

#Multivariate Mann-Kendall Test

require(trend)

data(hcb)

> plot(hcb)

 

 

s <- maxau[,"s"]

Q <- maxau[,"Q"]

cor.test(s,Q, meth="spearman")

 

partial.mk.test(s,Q)

partial.cor.trend.test(s,Q, "spearman")

 

#Cox and Stuart Trend test

cs.test(prectot)

 

s <- maxau[,"s"]

sens.slope(prectot)

 

dat=read.csv(file("clipboard"),header=T, sep="\t", row.names=1)

require(trend)

das=ts(dat, frequency=12,start=1961)

sea.sens.slope(das)

 

 

# Changement Dectection

#data(PagesData)

require(trend)

data=read.csv(file("clipboard"),header=T, sep="\t", row.names=1)

attach(data)

 

 

#Pettitt's test for single change-point detection

pettitt.test(dat)

 

 

#Buishand range test

require(trend)

(res <- br.test(prectot))

 

require(trend)

(res <- snh.test(Nile))

 

#Wallis and Moore Phase-Frequency test

frost <- ts(prectot, start=1961)

wm.test(frost)

#Bartels's test for randomness

bartels.test(prectot)

 

#Wald-Wolfowitz test for independence and stationarity

ww.test(prectot)

 

#Wallis and Moore Phase-Frequency test

wm.test(prectot)

 

#Buishand U test

#Nile

(res <- bu.test(prectot))

 

precto=ts(prectot, frequency=12, start=12)

par(mfrow=c(1,2))

plot(precto); plot(res)

 

 

#Standard Normal Homogeinity Test

require(trend)

#Nile

(res <- snh.test(prectot))

 

 

cor.test( x,y,alternative = c("two.sided", "less", "greater"),

use="complete.obs", method="kendall",conf.level = 0.95)

 

require(graphics)

x <- rnorm(50)

y <- runif(30)

# Do x and y come from the same distribution?

ks.test(x, y)

# Does x come from a shifted gamma distribution with shape 3 and rate 2?

ks.test(x+2, "pgamma", 3, 2) # two-sided, exact

ks.test(x+2, "pgamma", 3, 2, exact = FALSE)

ks.test(x+2, "pgamma", 3, 2, alternative = "gr")

 

#  The approach after Pettitt (1979)

#  is commonly applied to detect a single change-point

#  in hydrological series or climate series with continuous data.

#  It tests the H0: The T variables follow one or more

#  distributions that have the same location parameter (no change),

#  against the alternative: a change point exists.

 

The Pettitt-test is conducted in such a way

 

# Change-point Detection: Pettitt’s Test

library(package = dplyr, warn.conflicts = FALSE )

library(package = ggplot2, warn.conflicts = FALSE)

library(package = trend, warn.conflicts = FALSE)

 

# Load Data

data <- read.csv(file = 'c:/Users/pooya/Downloads/Torbat Heydariyeh - Daily.csv',

                 header = TRUE)

data=read.csv(file("clipboard"),header=T, sep="\t", row.names=1)

 

# Compactly Display the Structure of data

str(object = data)

 

 

# Return the First 10 Row of data

head(x = data, n = 10)

 

# Select Variable

Tave <- data %>%

  mutate(date = as.Date(x = paste0(Year, '-', Month, '-', Day), '%Y-%m-%d')) %>%

  select(date, t)

 

# Compactly Display the Structure of Tave

str(object = Tave)

# Return the First 10 Row of Tave

head(x = Tave, n = 10)

 

# Tave Summaries

summary(object = Tave)

 

# NA Remove

Tave <- na.omit(Tave)

# Plot

require(ggplot2)

ggplot(data = Tave, mapping = aes(x = date, y = prectot)) +

  geom_line()

 

# Pettitt’s test

pettittTest <- trend::pettitt.test(x = Tave[['prectot']])

print(pettittTest)

 

 

print(Tave[['date']][pettittTest$estimate])

 

 

# Plot

ggplot(data = Tave, mapping = aes(x = date, y = t)) +

  geom_line() +

  geom_vline(mapping = aes(xintercept = as.numeric(Tave[['date']][pettittTest$estimate])),

             linetype = 2,

             colour = "red",

             size = 2)

 

 

#  Break Point

#first generate some data that has monthly observations on 6 years of a seasonal

#process followed by 2 years of monthly that is not seasonal

 

#-------------------------------------------------------------------------------------

ry=read.table(file("clipboard"),header=T, sep="\t", row.names=1)

attach(ry)

str(ry)

 

prectot=ts(prectot, start=1961,end=2017)

plot(prectot, ylab="Precipitation (mm)", type="o",xlab="Year")

abline(h= mean(prectot),col='blue')

require(strucchange)

prectot=ts(prectot, start=1961,end=2017)

ocus.nile <- efp(prectot ~ 1, type = "OLS-CUSUM")

omus.nile <- efp(prectot ~ 1, type = "OLS-MOSUM")

rocus.nile <- efp(prectot ~ 1, type = "Rec-CUSUM")

 

 

opar <- par(mfrow=c(2,2))

plot(prectot, ylab="Precipitation (mm)", type="o",xlab="Year", main ="Annual Precipitation")

abline(h= mean(prectot),col='blue')

plot(ocus.nile, type="o",xlab="Year")

plot(omus.nile, type="o",xlab="Year")

plot(rocus.nile, type="o",xlab="Year")

require(strucchange)

attach(rm)

#first we are going to use the supF test to determine whether the parameters of #our monthly dummy variable regression are stable through the time series:

 

fs.days = Fstats(year~Dec+Feb+Mar+Apr+May+Jun+Jul+Aug+Sep+Oct+Nov,data=rm)

plot(fs.days)

 

#use the breakpoints functions to find BIC corresponding to optimal number of #breaks

bp.days = breakpoints(year~Dec+Feb+Mar+Apr+May,data=rm)

summary(bp.days)

x11()

plot(bp.days)

 

prectot=ts(prectot, start=1961, end=2017)

bp.nile <- breakpoints(prectot ~ 1)

nile.fac <- breakfactor(bp.nile, breaks = 1 )

fm1.nile <- lm(prectot~ nile.fac - 1)

plot(bp.nile)

 

prectot=ts(prectot, start=1961, end=2017)

opar <- par(mfrow=c(2,1), mar=c(2,2,0,2))

plot(ocus.nile, alt.boundary=F,main="")

abline(v= 1980, lty=2, col='red')

plot(prectot, ylab="")

abline(h= mean(prectot),col='blue')

abline(v= 1980, lty=2, col='red')

lines(ts(predict(fm1.nile),start=1980,freq=1), col='darkgreen',lwd=2)

par(opar)

 

 

 

plot(prectot)

 

## F statistics indicate one breakpoint

fs.nile <- Fstats(prectot ~ 1)

plot(fs.nile)

breakpoints(fs.nile)

lines(breakpoints(fs.nile))

 

## or

bp.nile <- breakpoints(prectot ~ 1)

summary(bp.nile)

 

## the BIC also chooses one breakpoint

plot(bp.nile)

breakpoints(bp.nile)

 

## confidence interval

ci.nile <- confint(bp.nile)

ci.nile

lines(ci.nile)

 

 

## fit null hypothesis model and model with 1 breakpoint

fm0 <- lm(prectot ~ 1)

fm1 <- lm(prectot ~ breakfactor(bp.nile, breaks = 1))

plot(prectot, xlab="Year", ylab="Preciptation (mm)")

lines(ts(fitted(fm0), start = 2006), col = 3)

lines(ts(fitted(fm1), start = 2006), col = 4)

lines(bp.nile)

abline(v=2006, lty=2, col=2)

 

## confidence intervals

ci.seat2 <- confint(bp.nile, breaks = 2)

ci.seat2

lines(ci.seat2)

 

 

#rm(list=ls())

data=read.csv(file("clipboard"),header=T, sep="\t")

str(data)

 

if(!require(mblm)){install.packages("mblm")}

if(!require(ggplot2)){install.packages("ggplot2")}

 

 

data$Year = round(data$Year)

data$Year = data$Year - 1961

 

 

library(mblm)

model= mblm(prectot ~ Year, data=data)

summary(model)

Sum = summary(model)$coefficients

 

library(ggplot2)

ggplot(data, aes(x=Year, y=prectot)) +

  geom_point() +

  geom_abline(intercept = Sum[1], slope = Sum[2], color="blue", size=1.2) +

  labs(x = "Years after 1961")

 

 

### Optional quantile regression follows

if(!require(quantreg)){install.packages("quantreg")}

library(quantreg)

 

model.q = rq(prectot ~ Year, data = data, tau = 0.5)

summary(model.q)

 

library(ggplot2)

model.null = rq(prectot ~ 1, data = data, tau = 0.5)

anova(model.q, model.null)

 

Sumq = summary(model.q)$coefficients

ggplot(data, aes(x=Year, y=prectot)) +

  geom_point() +

  geom_abline(intercept = Sumq[1], slope = Sumq[2], color="red", size=1.2) +

  labs(x = "Years after 1961")

 

 

1)Seasonality

###############

r=read.table(file("clipboard"),header=T, sep="\t", dec=".")

attach(r)

 

reg=lm(JFM~year)

summary(reg)

 

JFM=ts(JFM, start=1953, end=2017, freq=1)

plot(JFM, type="b", lty=3)

 

y = -0.25x + 519.54

p-value = .09054

 

 

rm=read.table(file("clipboard"),header=T, sep="\t", row.names=1)

getwd()

setwd("C:/Users/Lenovo/Documents")

write.table(rm, file="anclim.txt",row.names=F)

 

 

rm=read.table(file("clipboard"),header=T, sep="\t")

 

str(rm)

attach(rm)

 

 

#plot(0,0,xlim=c(1961,2017), ylim=c(0,300))

DJF=ts(DJF, start=1961,end=2017)

Y=ts(Y, start=1961, end=2017)

#par(new=TRUE)

#plot(DJF, type="o", xlab="", ylab="",axes=FALSE)

#par(new=TRUE)

#plot(Y, type="l", xlab="", ylab="", add=TRUE, axes=FALSE)

 

names(r)

reg=lm(SNO.tav~X,data=r)

summary(reg)

 

plot(X,MAM, type="o",lty=9, xlab="Year", ylab="Precipitation(mm)")

d=abline(reg)

legend("topright",c("precipitation(mm)","adjusted"),col=c("black","black"),cex=c(0.7,0.7),lty=c(9,1))

legend(locator(1),c("y= -0.323 x+677

P=.526

R²=.0073"),cex=0.65,bg="NA",box.col="white")

 

 

 

str(r)

library(climatol) # load the functions of the package

homogen("r", 1980, 2017, expl=TRUE)

#Missing Value

library(naniar)

vis_miss(r[,-1])

#dd2m("r", 1981, 2000, homog=TRUE)

#dahstat(’Ttest’, 1981, 2000, stat=’series’)

 

library(naniar)

require("UpSetR")

gg_miss_upset(r)

 

library(mice)

md.pattern(r)

 

library(reshape2)

library(ggplot2)

 

#Homogeneous

############

 

#Homogeneous

boxplot(data)

hist(data)

 

#median ajustement

pairs(datatable, col="blue", main="Scatterplots")

Y=cbind(LRY)

X=cbind(LRV, LRC, INT)

#

hist(Y, prob=TRUE, col = "blue", border = "black")

lines(density(Y))

#

OLSreg=lm(Y~X)

summary(OLSreg)

#

Qreg25=rq(Y~X, tau=0.25)

summary(Qreg25)

#

QR=rq(Y~X, tau=seq(0.2, 0.8, by=0.1))

sumQR=summary(QR)

#

anova(Qreg25, Qreg75)

#

QR=rq(Y~X, tau=seq(0.2, 0.8, by=0.1))

sumQR=summary(QR)

#

plot(sumQR)

 

 

 

# Check for Break point

 library("bfast")

 library("zoo")

 library("strucchange")

 

 plot(Nile, ylab="Annual Flow of the river Nile")

    abline(h= mean(Nile),col='blue')

 

 plot(merge(

 Nile = as.zoo(Nile),

 zoo(mean(Nile), time(Nile)),

 CUSUM = cumsum(Nile - mean(Nile)),

 zoo(0, time(Nile)),

 MOSUM = rollapply(Nile - mean(Nile), 15, sum),

 zoo(0, time(Nile))

 ), screen = c(1, 1, 2, 2, 3, 3), main = "", xlab = "Time",

 col = c(1, 4, 1, 4, 1, 4)

 ) 

 

plot(merge(

Nile = as.zoo(Nile),

zoo(c(NA, cumsum(head(Nile, -1))/1:99), time(Nile)),

CUSUM = cumsum(c(0, recresid(lm(Nile ~ 1)))),

zoo(0, time(Nile))

), screen = c(1, 1, 2, 2), main = "", xlab = "Time",

col = c(1, 4, 1, 4)

)

 

ocus.nile <- efp(Nile ~ 1, type = "OLS-CUSUM")

omus.nile <- efp(Nile ~ 1, type = "OLS-MOSUM")

rocus.nile <- efp(Nile ~ 1, type = "Rec-CUSUM")

 

 

opar <- par(mfrow=c(2,2))

plot(ocus.nile)

plot(omus.nile)

plot(rocus.nile)

par(opar)

 

 

plot(1870 + 5:95, sapply(5:95, function(i) {

before <- 1:i

after <- (i+1):100

res <- c(Nile[before] - mean(Nile[before]), Nile[after] - mean(Nile[after]))

sum(res^2)

}), type = "b", xlab = "Time", ylab = "RSS")

 

bp.nile <- breakpoints(Nile ~ 1)

nile.fac <- breakfactor(bp.nile, breaks = 1 )

fm1.nile <- lm(Nile ~ nile.fac - 1)

plot(bp.nile)

 

 

opar <- par(mfrow=c(2,1), mar=c(2,2,0,2))

plot(ocus.nile, alt.boundary=F,main="")

abline(v= 1898, lty=2, col='red')

plot(Nile, ylab="Annual Flow of the river Nile")

abline(h= mean(Nile),col='blue')

abline(v= 1898, lty=2, col='red')

lines(ts(predict(fm1.nile),start=1871,freq=1), col='darkgreen',lwd=2)

par(opar)

 

 

 

#TREND ANANLYSIS

#################

data=read.csv(file("clipboard"),header=T, sep="\t", row.names=1)

data

str(data)

 

#install.packages("trend")

require(trend)

 

#data(maxau)

prectot <- data[,"prectot"]

mk.test(prectot)

 

dat=read.csv(file("clipboard"),header=T, sep="\t", row.names=1)

data(nottem)

str(dat)

dat=round(dat,1)

dat1 <- ts(dat, frequency=12, start=c(1961,1))

plot.ts(dat1)

 

logsouvenirtimeseries <- log(dat1)

plot.ts(logsouvenirtimeseries)

 

require(trend)

smk.test(dat1)

 

#correlated witn precedent month

csmk.test(dat)

 

#Multivariate Mann-Kendall Test

require(trend)

data(hcb)

> plot(hcb)

 

 

s <- maxau[,"s"]

Q <- maxau[,"Q"]

cor.test(s,Q, meth="spearman")

 

partial.mk.test(s,Q)

partial.cor.trend.test(s,Q, "spearman")

 

#Cox and Stuart Trend test

cs.test(prectot)

 

s <- maxau[,"s"]

sens.slope(prectot)

 

dat=read.csv(file("clipboard"),header=T, sep="\t", row.names=1)

require(trend)

das=ts(dat, frequency=12,start=1961)

sea.sens.slope(das)

 

 

 

# Changement Dectection

#######################

#data(PagesData)

require(trend)

data=read.csv(file("clipboard"),header=T, sep="\t", row.names=1)

attach(data)

 

 

#Pettitt's test for single change-point detection

pettitt.test(dat)

 

 

#Buishand range test

require(trend)

(res <- br.test(prectot))

 

require(trend)

(res <- snh.test(Nile))

 

#Wallis and Moore Phase-Frequency test

frost <- ts(prectot, start=1961)

wm.test(frost)

 

 

#Bartels's test for randomness

bartels.test(prectot)

 

#Wald-Wolfowitz test for independence and stationarity

ww.test(prectot)

 

#Wallis and Moore Phase-Frequency test

wm.test(prectot)

 

 

#Buishand U test

#Nile

(res <- bu.test(prectot))

 

precto=ts(prectot, frequency=12, start=12)

par(mfrow=c(1,2))

plot(precto); plot(res)

 

 

#Standard Normal Homogeinity Test

require(trend)

#Nile

(res <- snh.test(prectot))

 

 

cor.test( x,y,alternative = c("two.sided", "less", "greater"),

use="complete.obs", method="kendall",conf.level = 0.95)

 

require(graphics)

x <- rnorm(50)

y <- runif(30)

# Do x and y come from the same distribution?

ks.test(x, y)

# Does x come from a shifted gamma distribution with shape 3 and rate 2?

ks.test(x+2, "pgamma", 3, 2) # two-sided, exact

ks.test(x+2, "pgamma", 3, 2, exact = FALSE)

ks.test(x+2, "pgamma", 3, 2, alternative = "gr")

 

 

#  The approach after Pettitt (1979)

#  is commonly applied to detect a single change-point

#  in hydrological series or climate series with continuous data.

#  It tests the H0: The T variables follow one or more

#  distributions that have the same location parameter (no change),

#  against the alternative: a change point exists.

 

The Pettitt-test is conducted in such a way

 

# Change-point Detection: Pettitt’s Test

library(package = dplyr, warn.conflicts = FALSE )

library(package = ggplot2, warn.conflicts = FALSE)

library(package = trend, warn.conflicts = FALSE)

 

# Load Data

data <- read.csv(file = 'c:/Users/pooya/Downloads/Torbat Heydariyeh - Daily.csv',

                 header = TRUE)

data=read.csv(file("clipboard"),header=T, sep="\t", row.names=1)

 

# Compactly Display the Structure of data

str(object = data)

 

 

# Return the First 10 Row of data

head(x = data, n = 10)

 

# Select Variable

Tave <- data %>%

  mutate(date = as.Date(x = paste0(Year, '-', Month, '-', Day), '%Y-%m-%d')) %>%

  select(date, t)

 

# Compactly Display the Structure of Tave

str(object = Tave)

 

# Return the First 10 Row of Tave

head(x = Tave, n = 10)

 

# Tave Summaries

summary(object = Tave)

 

# NA Remove

Tave <- na.omit(Tave)

 

# Plot

require(ggplot2)

ggplot(data = Tave, mapping = aes(x = date, y = prectot)) +

  geom_line()

 

 

# Pettitt’s test

pettittTest <- trend::pettitt.test(x = Tave[['prectot']])

print(pettittTest)

 

 

print(Tave[['date']][pettittTest$estimate])

 

 

# Plot

ggplot(data = Tave, mapping = aes(x = date, y = t)) +

  geom_line() +

  geom_vline(mapping = aes(xintercept = as.numeric(Tave[['date']][pettittTest$estimate])),

             linetype = 2,

             colour = "red",

             size = 2)

 

 

 

#  Break Point

#first generate some data that has monthly observations on 6 years of a seasonal

#process followed by 2 years of monthly that is not seasonal

 

#-------------------------------------------------------------------------------------

ry=read.table(file("clipboard"),header=T, sep="\t", row.names=1)

attach(ry)

str(ry)

 

prectot=ts(prectot, start=1961,end=2017)

plot(prectot, ylab="Precipitation (mm)", type="o",xlab="Year")

abline(h= mean(prectot),col='blue')

require(strucchange)

prectot=ts(prectot, start=1961,end=2017)

ocus.nile <- efp(prectot ~ 1, type = "OLS-CUSUM")

omus.nile <- efp(prectot ~ 1, type = "OLS-MOSUM")

rocus.nile <- efp(prectot ~ 1, type = "Rec-CUSUM")

 

 

opar <- par(mfrow=c(2,2))

plot(prectot, ylab="Precipitation (mm)", type="o",xlab="Year", main ="Annual Precipitation")

abline(h= mean(prectot),col='blue')

plot(ocus.nile, type="o",xlab="Year")

plot(omus.nile, type="o",xlab="Year")

plot(rocus.nile, type="o",xlab="Year")

 

 

 

require(strucchange)

attach(rm)

#first we are going to use the supF test to determine whether the parameters of #our monthly dummy variable regression are stable through the time series:

 

fs.days = Fstats(year~Dec+Feb+Mar+Apr+May+Jun+Jul+Aug+Sep+Oct+Nov,data=rm)

plot(fs.days)

 

#use the breakpoints functions to find BIC corresponding to optimal number of #breaks

bp.days = breakpoints(year~Dec+Feb+Mar+Apr+May,data=rm)

summary(bp.days)

x11()

plot(bp.days)

 

 

prectot=ts(prectot, start=1961, end=2017)

bp.nile <- breakpoints(prectot ~ 1)

nile.fac <- breakfactor(bp.nile, breaks = 1 )

fm1.nile <- lm(prectot~ nile.fac - 1)

plot(bp.nile)

 

prectot=ts(prectot, start=1961, end=2017)

opar <- par(mfrow=c(2,1), mar=c(2,2,0,2))

plot(ocus.nile, alt.boundary=F,main="")

abline(v= 1980, lty=2, col='red')

plot(prectot, ylab="")

abline(h= mean(prectot),col='blue')

abline(v= 1980, lty=2, col='red')

lines(ts(predict(fm1.nile),start=1980,freq=1), col='darkgreen',lwd=2)

par(opar)

 

 

plot(prectot)

 

## F statistics indicate one breakpoint

fs.nile <- Fstats(prectot ~ 1)

plot(fs.nile)

breakpoints(fs.nile)

lines(breakpoints(fs.nile))

 

## or

bp.nile <- breakpoints(prectot ~ 1)

summary(bp.nile)

 

## the BIC also chooses one breakpoint

plot(bp.nile)

breakpoints(bp.nile)

 

## confidence interval

ci.nile <- confint(bp.nile)

ci.nile

lines(ci.nile)

 

 

## fit null hypothesis model and model with 1 breakpoint

fm0 <- lm(prectot ~ 1)

fm1 <- lm(prectot ~ breakfactor(bp.nile, breaks = 1))

plot(prectot, xlab="Year", ylab="Preciptation (mm)")

lines(ts(fitted(fm0), start = 2006), col = 3)

lines(ts(fitted(fm1), start = 2006), col = 4)

lines(bp.nile)

abline(v=2006, lty=2, col=2)

 

## confidence intervals

ci.seat2 <- confint(bp.nile, breaks = 2)

ci.seat2

lines(ci.seat2)

 

 

#rm(list=ls())

data=read.csv(file("clipboard"),header=T, sep="\t")

str(data)

 

if(!require(mblm)){install.packages("mblm")}

if(!require(ggplot2)){install.packages("ggplot2")}

 

 

data$Year = round(data$Year)

data$Year = data$Year - 1961

 

 

library(mblm)

model= mblm(prectot ~ Year, data=data)

summary(model)

Sum = summary(model)$coefficients

 

library(ggplot2)

ggplot(data, aes(x=Year, y=prectot)) +

  geom_point() +

  geom_abline(intercept = Sum[1], slope = Sum[2], color="blue", size=1.2) +

  labs(x = "Years after 1961")

 

 

### Optional quantile regression follows

if(!require(quantreg)){install.packages("quantreg")}

library(quantreg)

 

model.q = rq(prectot ~ Year, data = data, tau = 0.5)

summary(model.q)

 

library(ggplot2)

model.null = rq(prectot ~ 1, data = data, tau = 0.5)

anova(model.q, model.null)

 

Sumq = summary(model.q)$coefficients

ggplot(data, aes(x=Year, y=prectot)) +

  geom_point() +

  geom_abline(intercept = Sumq[1], slope = Sumq[2], color="red", size=1.2) +

  labs(x = "Years after 1961")

 

#Calculus for Index CLIMATE CHANGE WITH R

#   0>   >>   <<  Importation Data

rm(list=ls())

djibclim=read.table(file("clipboard"),header=T, dec=".",sep="\t")

#File Change dir

#getwd()

#setwd("D:/Docs.Abdi-Basid.ADAN/14.Doc.LMME/A.ClimDjibouti/0.Data")

#save(djiclim, file = "djiclim.rda")

#load(file = "djibclim.rda")

#View(djibclim)

#str(djibclim);summary(djiclim)

#sapply(djibclim,sd)

#sapply(djibclim, mean, na.rm=TRUE)

#library(psych)

#describe(djiclim[,4:6])

#Different Plot Types:

#plot(r,type="l",col=2)

#plot(r,type="b")

#plot(r,type="p")

#plot(r,type="o")

#plot(r,type="h")

#plot(r,type="s", col=4)

#plot(r,type="c")

plot(r,type="o", main="", xlab="Year", ylab="Precipitation (mm)")

attach(r)

summary(lm(prectot~Year, data=r))

#

#   1>   >>   <<  source import Rclimtool

#File change dir

#getwd()

setwd("D:/Docs.Abdi-Basid.ADAN/14.Doc.LMME/2.Project 2")

source("RClimTool\\RClimTool.R")

#source("RClimTool\\RClimTool_v2.R")

#folder na gen, na, missing before homogeneous

#day month year  without label

#

#tmax = cargar(gfile(type = "open"))

#tmin = cargar(gfile(type = "open"))

#precip = cargar(gfile(type = "open"))

#

#plot(x = day, main = "Title", xlab = "Label for the x axis", ylab = "Label for the y axis")

#hist(x = day, main = "Title", xlab = "Label for the x axis", ylab = "Label for the y axis")

#boxplot(x = day, main = "Title", ylab = "Label for the y axis")

#

setwd("D:\\Docs.Abdi-Basid.ADAN\\14.Doc.LMME\\Rclimtrends\\R")

#library(climtrends)

source("climtrends.R")

names(djibclim)

attach(djibclim)

CraddockTest(Days[,Temp_Tn])

#

# Craddock test for testing the homogeneity of yearly average temperatures of Milan vs Turin

data(yearly.average.temperature.Turin.Milan)

craddock.result <- CraddockTest(yearly.average.temperature.Turin.Milan,2,3)

plot(yearly.average.temperature.Turin.Milan[,1],craddock.result,type='l',xlab='',ylab='Craddock')

#

# Example of inhomogeneity

tmp <- yearly.average.temperature.Turin.Milan

tmp[20:43,3] <- tmp[20:43,3]+1 # adding an error of +1 degree to Milan's yearly average temperatures from 1980

craddock.result <- CraddockTest(tmp,2,3)

dev.new()

plot(yearly.average.temperature.Turin.Milan[,1],craddock.result,type='l',xlab='',ylab='Craddock')

#

#Homogénéity Analysis

names(djiclim)

library(ggpubr)

ggqqplot(djiclim$Rain_Rrr)

#

library("car")

qqPlot(djiclim$Rain_Rrr)

#

#shapiro.test(rnorm(100, mean = 5, sd = 3))

shapiro.test(runif(100, min = 2, max = 4))

#require(stats)

shapiro.test(djibclim$Rain_Rrr)

#

library(tseries)

jarque.bera.test(djibclim$Rain_Rrr)

#

#Kolmogorov-Smirnov Tests

t.test(djibclim$Rain_Rrr)

wilcox.test(djibclim$Rain_Rrr)

ks.test(djibclim$Rain_Rrr)

#

#

#

#   2>   >>   << Program Rclim Dex

#library("RClimDex")

#rclimdex.start()

getwd()

setwd("D:/Docs.Abdi-Basid.ADAN/14.Doc.LMME/2.Project 2/Climat/RclimDex")

source("1.0.RclimDex.R.txt")

require("climdex.pcic")

require("Kendall")

require(tcltk)

#year month day without label

#I+1  i-1

#

#   3>   >>   << Program ClimPact2

#Download the above zip file to your computer

#(https://github.com/ARCCSS-extremes/climpact2/archive/master.zip)

# and extract it. This will create a directory named climpact2-master.

#In Windows: open R and select "File -> Change dir..."

#and select the climpact2-master directory created in step 1.

#Then type source('climpact2.GUI.r').

#(NOTE: If nothing happens, try run the additional command startss()).

#In Linux/macOS: change to the climpact2-master directory created in step 1,

#then open R in a terminal window and type source('climpact2.GUI.r').

#The first time Climpact2 is run it will install required R packages.

#This will likely require you to select a mirror from which to download.

getwd()

setwd("D:\\Docs.Abdi-Basid.ADAN\\14.Doc.LMME\\2.Project 2\\Climat\\RClimPact\\climpact2-master")

source("climpact2.GUI.r")

#

#

#   4>   >>   << Program RHtest

#setwd("D:\\Docs.Abdi-Basid.ADAN\\Doc.LMME\\RHtest")

#source ("RHtest_V4.r") 

#source ("RHtest.txt") 

setwd("D:\\Docs.Abdi-Basid.ADAN\\14.Doc.LMME\\RHtest\\RHtests-master")

source("RHtestsV4_20180301.r")

StartGUI()

#   4'>   >>   << Program RHtest

library(PCICt)

## Create a climdexInput object from some data already loaded in and

## ready to go.

## Parse the dates into PCICt.

tmax.dates <- as.PCICt(do.call(paste, ec.1018935.tmax[,c("year",

"jday")]), format="%Y %j", cal="gregorian")

tmin.dates <- as.PCICt(do.call(paste, ec.1018935.tmin[,c("year",

"jday")]), format="%Y %j", cal="gregorian")

prec.dates <- as.PCICt(do.call(paste, ec.1018935.prec[,c("year",

"jday")]), format="%Y %j", cal="gregorian")

## Load the data in.

ci <- climdexInput.raw(ec.1018935.tmax$MAX_TEMP,

ec.1018935.tmin$MIN_TEMP, ec.1018935.prec$ONE_DAY_PRECIPITATION,

tmax.dates, tmin.dates, prec.dates, base.range=c(1971, 2000))

## Create an annual timeseries of the CDD index.

cdd <- climdex.cdd(ci)

#

## Create a climdexInput object from some data already loaded in and

## ready to go.

## Parse the dates into PCICt.

tmax.dates <- as.PCICt(do.call(paste, ec.1018935.tmax[,c("year",

"jday")]), format="%Y %j", cal="gregorian")

tmin.dates <- as.PCICt(do.call(paste, ec.1018935.tmin[,c("year",

"jday")]), format="%Y %j", cal="gregorian")

prec.dates <- as.PCICt(do.call(paste, ec.1018935.prec[,c("year",

"jday")]), format="%Y %j", cal="gregorian")

## Load the data in.

ci <- climdexInput.raw(ec.1018935.tmax$MAX_TEMP,

ec.1018935.tmin$MIN_TEMP, ec.1018935.prec$ONE_DAY_PRECIPITATION,

tmax.dates, tmin.dates, prec.dates, base.range=c(1971, 2000))

## Create an annual timeseries of the cold spell duration index.

csdi <- climdex.csdi(ci)

#

#

#

## Create a climdexInput object from some data already loaded in and

## ready to go.

## Parse the dates into PCICt.

tmax.dates <- as.PCICt(do.call(paste, ec.1018935.tmax[,c("year",

"jday")]), format="%Y %j", cal="gregorian")

tmin.dates <- as.PCICt(do.call(paste, ec.1018935.tmin[,c("year",

"jday")]), format="%Y %j", cal="gregorian")

prec.dates <- as.PCICt(do.call(paste, ec.1018935.prec[,c("year",

"jday")]), format="%Y %j", cal="gregorian")

## Load the data in.

ci <- climdexInput.raw(ec.1018935.tmax$MAX_TEMP,

ec.1018935.tmin$MIN_TEMP, ec.1018935.prec$ONE_DAY_PRECIPITATION,

tmax.dates, tmin.dates, prec.dates, base.range=c(1971, 2000))

## Create an annual timeseries of the CWD index.

cdd <- climdex.cdd(ci)

 

#

#   5>   >>   <<  ANCLIM

 

at=read.table(file("clipboard"),header=T, sep="\t")

write(at, file = "data.txt", sep = " ")

 

#   5>   >>   <<  Geoclim

 

 

#   6>   >>   <<  GeoDa

 

 

 

 

 

 

 

 

 

 

 

 

 

# ANNEXE

########

 

require(ggplot2)

ggplot(data = quakes) +

       geom_density(mapping = aes(x = mag, fill = region), alpha = 0.5)

 

library("RColorBrewer")

display.brewer.all()

 

 

 

# Missing Data Imputation

#########################

library(mice)

Mousey = mice(Greenwich)

Greenwich = complete(Mousey)

 

misp=read.csv(file("clipboard"), header=T, sep="\t")

str(misp)

misp=misp[,-c(1:3)]

library(naniar)

vis_miss(misp)

 

 

#Prediction

###########

 

#Import Data

data=read.table(file("clipboard"), header=T, sep="\t", dec=".")

str(data)

summary(data)

library(forecast)

 

fit.series3 <- auto.arima(prectot)

fcast.series3 <- forecast(fit.series3)

plot(fcast.series3)

 

ts_passengers = ts(prectot,start=1961,end=2017,frequency=1)

plot(ts_passengers,xlab="", ylab="")

 

 

m_ets = ets(ts_passengers)

f_ets = forecast(m_ets, h=24) # forecast 24 months into the future

plot(f_ets, type="o", xlab="Year", ylab="Precipitaion(mm)")

 

 

m_tbats = tbats(ts_passengers)

f_tbats = forecast(m_tbats, h=24)

plot(f_tbats)

 

 

m_aa = auto.arima(ts_passengers)

f_aa = forecast(m_aa, h=24)

plot(f_aa)

 

 

barplot(c(ETS=m_ets$aic, ARIMA=m_aa$aic, TBATS=m_tbats$AIC),

    col="light blue",

    ylab="AIC")

 

 

#https://a-little-book-of-r-for-time-series.readthedocs.io/en/latest/src/timeseries.html

#https://otexts.com/fpp2/arima-r.html

 

 

prectot %>% stl(s.window='periodic') %>% seasadj() -> prectot

autoplot(prectot)

 

#Steps to be followed for ARIMA modeling:

#1. Exploratory analysis

#2. Fit the model

#3. Diagnostic measures

tsData = ts(prectot, start = c(1961,1), end=c(2017,1),frequency = 12)

 

#1. Exploratory analysis

components.ts = decompose(tsData)

plot(components.ts)

 

library("fUnitRoots")

urkpssTest(tsData, type = c("tau"), lags = c("short"),use.lag = NULL, doplot = TRUE)

tsstationary = diff(tsData, differences=1)

plot(tsstationary)

 

acf(tsstationary, lag.max=34)

pacf(tsstationary, lag.max=34)

 

fitARIMA <- arima(tsData, order=c(1,1,1),seasonal = list(order = c(1,0,0), period = 12),method="ML")

library(lmtest)

coeftest(fitARIMA)

confint(fitARIMA)

 

acf(fitARIMA$residuals)

library(FitAR)

LjungBoxTest (fitARIMA$residuals,k=2,StartLag=1)

 

qqnorm(fitARIMA$residuals)

qqline(fitARIMA$residuals)

auto.arima(tsData, trace=TRUE)

 

predict(fitARIMA,n.ahead = 20)

require(forecast)

futurVal <- forecast.Arima(fitARIMA,h=10, level=c(99.5))

plot.forecast(futurVal)

plot(fitARIMA,h=10, level=c(99.5))

 

 

A=auto.arima(prectot)

autoplot(forecast(A))

 

 

autoplot(forecast(prectot), n=100)

 

library(forecast)

AutoArimaModel=auto.arima(prectot)

AutoArimaModel

 

plot(AutoArimaModel)

require(tseries); require(astsa)

acf(prectot)

pacf(prectot)

 

precto=ts(prectot, frequency=12, start=1961, end=2017)

 

fit = arima(precto, order = c(0, 0, 0))

pred = predict(fit, n.ahead = 50)

 

 

fit = arima(log(precto), c(1, 1, 0), seasonal = list(order = c(0, 1, 1), period = 12))

pred <- predict(fit, n.ahead = 10*12)

ts.plot(precto,exp(pred$pred), log = "y", lty = c(1,3))

 

od <- options(digits = 5) # avoid too much spurious accuracy

q=predict(arima(lh, order = c(3,0,0)), n.ahead = 12)

 

 

(fit <- arima(USAccDeaths, order = c(0,1,1),

              seasonal = list(order = c(0,1,1))))

 

predict(fit, n.ahead = 6)

 

 

 

ts(vector, start=, end=, frequency=)

# subset the time series (June 2014 to December 2014)

myts2 <- window(myts, start=c(2014, 6), end=c(2014, 12))

 

# plot series

plot(myts)

 

# Seasonal decomposition

fit <- stl(myts, s.window="period")

plot(fit)

 

# additional plots

monthplot(myts)

library(forecast)

seasonplot(myts)

 

 

# simple exponential - models level

fit <- HoltWinters(myts, beta=FALSE, gamma=FALSE)

# double exponential - models level and trend

fit <- HoltWinters(myts, gamma=FALSE)

# triple exponential - models level, trend, and seasonal components

fit <- HoltWinters(myts)

 

# predictive accuracy

library(forecast)

accuracy(fit)

 

# predict next three future values

library(forecast)

forecast(fit, 3)

plot(forecast(fit, 3))

 

 

# fit an ARIMA model of order P, D, Q

fit <- arima(myts, order=c(p, d, q)

 

# predictive accuracy

library(forecast)

accuracy(fit)

 

# predict next 5 observations

library(forecast)

forecast(fit, 5)

plot(forecast(fit, 5))

 

 

library(forecast)

# Automated forecasting using an exponential model

fit <- ets(myts)

 

# Automated forecasting using an ARIMA model

fit <- auto.arima(myts)

 

 

#The forecast package provides functions for the automatic selection

#of exponential and ARIMA models. The ets() function

#supports both additive and multiplicative models.

#The auto.arima() function can handle both seasonal

#and nonseasonal ARIMA models. Models are chosen to maximize one of

#several fit criteria.

 

#lag(ts, k)       lagged version of time series, shifted back k observations

#diff(ts, differences=d)          difference the time series d times

#ndiffs(ts)      Number of differences required to achieve stationarity (from the forecast package)

#acf(ts)           autocorrelation function

#pacf(ts)         partial autocorrelation function

#adf.test(ts)   Augemented Dickey-Fuller test. Rejecting the null hypothesis suggests that a time series is stationary (from the tseries package)

#Box.test(x, type="Ljung-Box")        Pormanteau test that observations in vector or time series x are independent

 

 

x <- c(32,64,96,118,126,144,152.5,158) 

y <- c(99.5,104.8,108.5,100,86,64,35.3,15)

 

 

x <- 1:10

y <- x + c(-0.5,0.5)

 

plot(x,y, xlim=c(0,11), ylim=c(-1,12))

 

fit1 <- lm( y~offset(x) -1 )

fit2 <- lm( y~x )

fit3 <- lm( y~poly(x,3) )

fit4 <- lm( y~poly(x,9) )

library(splines)

fit5 <- lm( y~ns(x, 3) )

fit6 <- lm( y~ns(x, 9) )

 

fit7 <- lm( y ~ x + cos(x*pi) )

 

xx <- seq(0,11, length.out=250)

lines(xx, predict(fit1, data.frame(x=xx)), col='blue')

lines(xx, predict(fit2, data.frame(x=xx)), col='green')

lines(xx, predict(fit3, data.frame(x=xx)), col='red')

lines(xx, predict(fit4, data.frame(x=xx)), col='purple')

lines(xx, predict(fit5, data.frame(x=xx)), col='orange')

lines(xx, predict(fit6, data.frame(x=xx)), col='grey')

lines(xx, predict(fit7, data.frame(x=xx)), col='black')

 

 

polyfit <- function(i) x <- AIC(lm(y~poly(x,i)))

as.integer(optimize(polyfit,interval = c(1,length(x)-1))$minimum)

 

for (i in 2:length(x)-1) print(polyfit(i))

 

 

lm(y ~ x + I(x^2) + I(x^3))

lm(y ~ poly(x, 3, raw=TRUE))

 

 

 

#Reseau Neurone Prediction

 

# reseaux neurone  application sous R

 

library(MASS)

library(nnet)

 

# apprentissage

nnet.reg=nnet(O3obs~.,data=datappr,size=5,decay=1,linout=TRUE,maxit=500)

summary(nnet.reg)

 

 

library(e1071)

plot(tune.nnet(O3obs~.,data=datappr,size=c(2,3,4),

decay=c(1,2,3),maxit=200,linout=TRUE))

plot(tune.nnet(O3obs~.,data=datappr,size=4:5,decay=1:10)

 

nnet.reg=nnet(O3obs~.,data=datappr,size=3,decay=2,

linout=TRUE,maxit=200)

# calcul et graphe des résidus

fit.nnetr=predict(nnet.reg,data=datappr)

res.nnetr=fit.nnetr-datappr[,"O3obs"]

plot.res(fit.nnetr,res.nnetr)

 

 

# apprentissage

nnet.dis=nnet(DepSeuil~.,data=datappq,size=5,decay=1)

summary(nnet.reg)

#matrice de confusion

table(nnet.dis$fitted.values>0.5,datappq$DepSeuil)

 

 

function(formula, data, size, niter = 1,

nplis = 10, decay = 0, maxit = 100)

{

n = nrow(data)

tmc=0

un = rep(1, n)

ri = sample(nplis, n, replace = T)

cat(" k= ")

for(i in sort(unique(ri))) {

cat(" ", i, sep = "")

for(rep in 1:niter) {

learn = nnet(formula, data[ri != i, ],

size = size,trace = F, decay = decay,

maxit = maxit)

tmc = tmc + sum(un[(data$DepSeuil[ri == i]

== "TRUE") != (predict(learn, data[ri == i,

]) > 0.5)])

}

}

cat("\n", "Taux de mal classes")

tmc/(niter * length(unique(ri)) * n)

}

 

CVnn(DepSeuil~.,data=datappq,size=7, decay=0)

...

# exécuter pour différentes valeur du decay

 

pred.nnetr=predict(nnet.reg,newdata=datestr)

pred.nnetq=predict(nnet.dis,newdata=datestq)

# Erreur quadratique moyenne de prévision

sum((pred.nnetr-datestr[,"O3obs"])^2)/nrow(datestr)

# Matrice de confusion pour la prévision du

# dépassement de seuil (régression)

table(pred.nnetr>150,datestr[,"O3obs"]>150)

# Même chose pour la discrimination

table(pred.nnetq>0.5,datestq[,"DepSeuil"])

 

library(ROCR)

rocnnetr=pred.nnetr/300

prednnetr=prediction(rocnnetr,datestq$DepSeuil)

perfnnetr=performance(prednnetr,"tpr","fpr")

rocnnetq=pred.nnetq

prednnetq=prediction(rocnnetq,datestq$DepSeuil)

perfnnetq=performance(prednnetq,"tpr","fpr")

# tracer les courbes ROC en les superposant

# pour mieux comparer

plot(perflogit,col=1)

plot(perfnnetr,col=2,add=TRUE)

plot(perfnnetq,col=3,add=TRUE)

 

 

 

 

#Extreme Value Analysis for Climate Change

r=read.table(file("clipboard"),header=TRUE,sep="\t", row.names=1)

str(r)

attach(r)

prcp=precipitation

moy=mean(precipitation)

 

sigma=sqrt(1/(length(precipitation)-1))*sum((precipitation-moy)^2)

k=(sigma/moy)^-1.086

 

 

x=precipitation

GEV=(1/sigma)*exp(- (1+(k*(x-moy)/sigma))^-(1/k)  )*(1+(k*(x-moy)/sigma ))^(-1-(1/k))

GEV

barplot(GEV,type="o")

GP=(1/sigma)*(1+(k*((x-moy)/sigma)) ^(-1-(1/k)))

GP

plot(GP,type="o")

GEV/GP

 

 

 

 

# load packages

library(extRemes)

library(xts)

 

# get data from eHYD

ehyd_url <- "http://ehyd.gv.at/eHYD/MessstellenExtraData/nlv?id=105700&file=2"

precipitation_xts <- read_ehyd(ehyd_url)

 

# mean residual life plot:

mrlplot(precipitation_xts, main="Mean Residual Life Plot")

# The mean residual life plot depicts the Thresholds (u) vs Mean Excess flow.

# The idea is to ?nd the lowest threshold where the plot is nearly linear;

# taking into account the 95% con?dence bounds.

 

# fitting the GPD model over a range of thresholds

threshrange.plot(precipitation_xts, r = c(30, 45), nint = 16)

# ismev implementation is faster:

# ismev::gpd.fitrange(precipitation_xts, umin=30, umax=45, nint = 16)

# set threshold

th <- 40

 

# maximum likelihood estimation

pot_mle <- fevd(as.vector(precipitation_xts), method = "MLE", type="GP", threshold=th)

# diagnostic plots

plot(pot_mle)

rl_mle <- return.level(pot_mle, conf = 0.05, return.period= c(2,5,10,20,50,100))

 

# L-moments estimation

pot_lmom <- fevd(as.vector(precipitation_xts), method = "Lmoments", type="GP", threshold=th)

# diagnostic plots

plot(pot_lmom)

rl_lmom <- return.level(pot_lmom, conf = 0.05, return.period= c(2,5,10,20,50,100))

 

# return level plots

par(mfcol=c(1,2))

# return level plot w/ MLE

plot(pot_mle, type="rl",

     main="Return Level Plot for Oberwang w/ MLE",

     ylim=c(0,200), pch=16)

loc <- as.numeric(return.level(pot_mle, conf = 0.05,return.period=100))

segments(100, 0, 100, loc, col= 'midnightblue',lty=6)

segments(0.01,loc,100, loc, col='midnightblue', lty=6)

 

# return level plot w/ LMOM

plot(pot_lmom, type="rl",

     main="Return Level Plot for Oberwang w/ L-Moments",

     ylim=c(0,200))

loc <- as.numeric(return.level(pot_lmom, conf = 0.05,return.period=100))

segments(100, 0, 100, loc, col= 'midnightblue',lty=6)

segments(0.01,loc,100, loc, col='midnightblue', lty=6)

 

# comparison of return levels

results <- t(data.frame(mle=as.numeric(rl_mle),

                        lmom=as.numeric(rl_lmom)))

colnames(results) <- c(2,5,10,20,50,100)

round(results,1)

 

 

 

 

 

 

library(extRemes)

library(ismev)

 

# Southwest-England rainfall data from Coles

data(rain)

 

# simple threshold (usually one should not use fixed quantiles)

u < - quantile(rain, 0.9)

 

# fit GP using Bayesian parameter estimation

x <- fevd(rain, threshold = u , type='GP', method='Bayesian')

 

 

 

 

B1.fit=fevd(B1[,2], data=B1, threshold=B1.th, type="GEV",method ="MLE", units="KVA",time.units="seconds",period = "Seconds")

B1.fit1=fevd(B1[,2], data=B1, threshold=B1.th,type="GP",method ="MLE", units="KVA",time.units="seconds",period = "Seconds")

B1.fit2=fevd(B1[,2], data=B1, threshold=B1.th,type="Gumbel", method ="MLE",units="KVA",time.units="seconds",period = "Seconds")

B1.fit3=fevd(B1[,2], data=B1, threshold=B1.th,type="Exponential",method ="MLE", units="KVA",time.units="seconds",period = "Seconds")

 

fit.AIC=summary(B1.fit, silent=TRUE)$AIC

fit1.AIC=summary(B1.fit1, silent=TRUE)$AIC

fit2.AIC=summary(B1.fit2, silent=TRUE)$AIC

fit3.AIC=summary(B1.fit3, silent=TRUE)$AIC

 

fit.AIC

# [1] 39976258

fit1.AIC

# [1] 466351.5

fit2.AIC

# [1] 13934878

fit3.AIC

# [1] 466330.8

 

plot(B1.fit)

plot(B1.fit1)

plot(B1.fit2)

plot(B1.fit3)

 

 

 

 

 

 

 

 

 

 

rm(list=ls())

rm(list=ls())

rm(list=ls())

 

 

rm=read.table(file("clipboard"),header=T, sep="\t", row.names=1)

getwd()

setwd("C:/Users/Lenovo/Documents")

write.table(rm, file="anclim.txt",row.names=F)

 

 

rm=read.table(file("clipboard"),header=T, sep="\t")

 

str(rm)

attach(rm)

 

 

#plot(0,0,xlim=c(1961,2017), ylim=c(0,300))

DJF=ts(DJF, start=1961,end=2017)

Y=ts(Y, start=1961, end=2017)

#par(new=TRUE)

#plot(DJF, type="o", xlab="", ylab="",axes=FALSE)

#par(new=TRUE)

#plot(Y, type="l", xlab="", ylab="", add=TRUE, axes=FALSE)

 

names(r)

reg=lm(SNO.tav~X,data=r)

summary(reg)

 

plot(X,MAM, type="o",lty=9, xlab="Year", ylab="Precipitation(mm)")

d=abline(reg)

legend("topright",c("precipitation(mm)","adjusted"),col=c("black","black"),cex=c(0.7,0.7),lty=c(9,1))

legend(locator(1),c("y= -0.323 x+677

P=.526

R²=.0073"),cex=0.65,bg="NA",box.col="white")

 

 

 

str(r)

library(climatol) # load the functions of the package

homogen("r", 1980, 2017, expl=TRUE)

#Missing Value

library(naniar)

vis_miss(r[,-1])

#dd2m("r", 1981, 2000, homog=TRUE)

#dahstat(’Ttest’, 1981, 2000, stat=’series’)

 

library(naniar)

require("UpSetR")

gg_miss_upset(r)

 

library(mice)

md.pattern(r)

library(reshape2)

library(ggplot2

Abdi-Basid Analysis

abdi-basid@outlook.com

 

 

#Import Data

data=read.table(file("clipboard"), header=T, sep="\t", dec=".")

str(data)

summary(data)

#Homogeneous

boxplot(data)

hist(data)

 

#median ajustement

pairs(datatable, col="blue", main="Scatterplots")

Y=cbind(LRY)

X=cbind(LRV, LRC, INT)

#

hist(Y, prob=TRUE, col = "blue", border = "black")

lines(density(Y))

#

OLSreg=lm(Y~X)

summary(OLSreg)

#

Qreg25=rq(Y~X, tau=0.25)

summary(Qreg25)

#

QR=rq(Y~X, tau=seq(0.2, 0.8, by=0.1))

sumQR=summary(QR)

#

anova(Qreg25, Qreg75)

#

QR=rq(Y~X, tau=seq(0.2, 0.8, by=0.1))

sumQR=summary(QR)

#

plot(sumQR)

 

# Check for Break point

 library("bfast")

 library("zoo")

 library("strucchange")

 

ts(vector, start=, end=, frequency=)

# subset the time series (June 2014 to December 2014)

myts2 <- window(myts, start=c(2014, 6), end=c(2014, 12))

 

# plot series

plot(myts)

 

# Seasonal decomposition

fit <- stl(myts, s.window="period")

plot(fit)

 

# additional plots

monthplot(myts)

library(forecast)

seasonplot(myts)

 

 

# simple exponential - models level

fit <- HoltWinters(myts, beta=FALSE, gamma=FALSE)

# double exponential - models level and trend

fit <- HoltWinters(myts, gamma=FALSE)

# triple exponential - models level, trend, and seasonal components

fit <- HoltWinters(myts)

 

# predictive accuracy

library(forecast)

accuracy(fit)

 

# predict next three future values

library(forecast)

forecast(fit, 3)

plot(forecast(fit, 3))

 

 

# fit an ARIMA model of order P, D, Q

fit <- arima(myts, order=c(p, d, q)

 

# predictive accuracy

library(forecast)

accuracy(fit)

 

# predict next 5 observations

library(forecast)

forecast(fit, 5)

plot(forecast(fit, 5))

 

 

library(forecast)

# Automated forecasting using an exponential model

fit <- ets(myts)

 

# Automated forecasting using an ARIMA model

fit <- auto.arima(myts)

 

 

#The forecast package provides functions for the automatic selection

#of exponential and ARIMA models. The ets() function

#supports both additive and multiplicative models.

#The auto.arima() function can handle both seasonal

#and nonseasonal ARIMA models. Models are chosen to maximize one of

#several fit criteria.

 

#lag(ts, k)       lagged version of time series, shifted back k observations

#diff(ts, differences=d)          difference the time series d times

#ndiffs(ts)      Number of differences required to achieve stationarity (from the forecast package)

#acf(ts)           autocorrelation function

#pacf(ts)         partial autocorrelation function

#adf.test(ts)   Augemented Dickey-Fuller test. Rejecting the null hypothesis suggests that a time series is stationary (from the tseries package)

#Box.test(x, type="Ljung-Box")        Pormanteau test that observations in vector or time series x are independent

 

 plot(Nile, ylab="Annual Flow of the river Nile")

    abline(h= mean(Nile),col='blue')

 

 plot(merge(

 Nile = as.zoo(Nile),

 zoo(mean(Nile), time(Nile)),

 CUSUM = cumsum(Nile - mean(Nile)),

 zoo(0, time(Nile)),

 MOSUM = rollapply(Nile - mean(Nile), 15, sum),

 zoo(0, time(Nile))

 ), screen = c(1, 1, 2, 2, 3, 3), main = "", xlab = "Time",

 col = c(1, 4, 1, 4, 1, 4)

 ) 

 

plot(merge(

Nile = as.zoo(Nile),

zoo(c(NA, cumsum(head(Nile, -1))/1:99), time(Nile)),

CUSUM = cumsum(c(0, recresid(lm(Nile ~ 1)))),

zoo(0, time(Nile))

), screen = c(1, 1, 2, 2), main = "", xlab = "Time",

col = c(1, 4, 1, 4)

)

 

ocus.nile <- efp(Nile ~ 1, type = "OLS-CUSUM")

omus.nile <- efp(Nile ~ 1, type = "OLS-MOSUM")

rocus.nile <- efp(Nile ~ 1, type = "Rec-CUSUM")

 

 

opar <- par(mfrow=c(2,2))

plot(ocus.nile)

plot(omus.nile)

plot(rocus.nile)

par(opar)

 

plot(1870 + 5:95, sapply(5:95, function(i) {

before <- 1:i

after <- (i+1):100

res <- c(Nile[before] - mean(Nile[before]), Nile[after] - mean(Nile[after]))

sum(res^2)

}), type = "b", xlab = "Time", ylab = "RSS")

 

 

bp.nile <- breakpoints(Nile ~ 1)

nile.fac <- breakfactor(bp.nile, breaks = 1 )

fm1.nile <- lm(Nile ~ nile.fac - 1)

plot(bp.nile)

 

 

 

 

opar <- par(mfrow=c(2,1), mar=c(2,2,0,2))

plot(ocus.nile, alt.boundary=F,main="")

abline(v= 1898, lty=2, col='red')

plot(Nile, ylab="Annual Flow of the river Nile")

abline(h= mean(Nile),col='blue')

abline(v= 1898, lty=2, col='red')

lines(ts(predict(fm1.nile),start=1871,freq=1), col='darkgreen',lwd=2)

par(opar)

 

 

x <- c(32,64,96,118,126,144,152.5,158) 

y <- c(99.5,104.8,108.5,100,86,64,35.3,15)

 

 

x <- 1:10

y <- x + c(-0.5,0.5)

 

plot(x,y, xlim=c(0,11), ylim=c(-1,12))

 

fit1 <- lm( y~offset(x) -1 )

fit2 <- lm( y~x )

fit3 <- lm( y~poly(x,3) )

fit4 <- lm( y~poly(x,9) )

library(splines)

fit5 <- lm( y~ns(x, 3) )

fit6 <- lm( y~ns(x, 9) )

 

fit7 <- lm( y ~ x + cos(x*pi) )

 

xx <- seq(0,11, length.out=250)

lines(xx, predict(fit1, data.frame(x=xx)), col='blue')

lines(xx, predict(fit2, data.frame(x=xx)), col='green')

lines(xx, predict(fit3, data.frame(x=xx)), col='red')

lines(xx, predict(fit4, data.frame(x=xx)), col='purple')

lines(xx, predict(fit5, data.frame(x=xx)), col='orange')

lines(xx, predict(fit6, data.frame(x=xx)), col='grey')

lines(xx, predict(fit7, data.frame(x=xx)), col='black')

 

 

polyfit <- function(i) x <- AIC(lm(y~poly(x,i)))

as.integer(optimize(polyfit,interval = c(1,length(x)-1))$minimum)

 

for (i in 2:length(x)-1) print(polyfit(i))

 

 

lm(y ~ x + I(x^2) + I(x^3))

lm(y ~ poly(x, 3, raw=TRUE))

 

 

 

 

 

 

 

 

 

 

library(mice)

Mousey = mice(Greenwich)

Greenwich = complete(Mousey)

 

misp=read.csv(file("clipboard"), header=T, sep="\t")

str(misp)

misp=misp[,-c(1:3)]

library(naniar)

vis_miss(misp)

 

 

#                       PROGRAM CLIMATE CHANGE WITH R

#                                                   abdi-basid@outlook.com

#

#

#

#   0>   >>   <<  Importation Data

rm(list=ls())

djibclim=read.table(file("clipboard"),header=T, dec=".",sep="\t")

#File Change dir

#getwd()

#setwd("D:/Docs.Abdi-Basid.ADAN/14.Doc.LMME/A.ClimDjibouti/0.Data")

#save(djiclim, file = "djiclim.rda")

#load(file = "djibclim.rda")

#View(djibclim)

#str(djibclim);summary(djiclim)

#sapply(djibclim,sd)

#sapply(djibclim, mean, na.rm=TRUE)

#library(psych)

#describe(djiclim[,4:6])

#Different Plot Types:

#plot(r,type="l",col=2)

#plot(r,type="b")

#plot(r,type="p")

#plot(r,type="o")

#plot(r,type="h")

#plot(r,type="s", col=4)

#plot(r,type="c")

plot(r,type="o", main="", xlab="Year", ylab="Precipitation (mm)")

attach(r)

summary(lm(prectot~Year, data=r))

#

#

#   1>   >>   <<  source import Rclimtool

#File change dir

#getwd()

setwd("D:/Docs.Abdi-Basid.ADAN/14.Doc.LMME/2.Project 2")

source("RClimTool\\RClimTool.R")

#source("RClimTool\\RClimTool_v2.R")

#folder na gen, na, missing before homogeneous

#day month year  without label

#

#tmax = cargar(gfile(type = "open"))

#tmin = cargar(gfile(type = "open"))

#precip = cargar(gfile(type = "open"))

#

#plot(x = day, main = "Title", xlab = "Label for the x axis", ylab = "Label for the y axis")

#hist(x = day, main = "Title", xlab = "Label for the x axis", ylab = "Label for the y axis")

#boxplot(x = day, main = "Title", ylab = "Label for the y axis")

#

setwd("D:\\Docs.Abdi-Basid.ADAN\\14.Doc.LMME\\Rclimtrends\\R")

#library(climtrends)

source("climtrends.R")

names(djibclim)

attach(djibclim)

CraddockTest(Days[,Temp_Tn])

#

# Craddock test for testing the homogeneity of yearly average temperatures of Milan vs Turin

data(yearly.average.temperature.Turin.Milan)

craddock.result <- CraddockTest(yearly.average.temperature.Turin.Milan,2,3)

plot(yearly.average.temperature.Turin.Milan[,1],craddock.result,type='l',xlab='',ylab='Craddock')

#

# Example of inhomogeneity

tmp <- yearly.average.temperature.Turin.Milan

tmp[20:43,3] <- tmp[20:43,3]+1 # adding an error of +1 degree to Milan's yearly average temperatures from 1980

craddock.result <- CraddockTest(tmp,2,3)

dev.new()

plot(yearly.average.temperature.Turin.Milan[,1],craddock.result,type='l',xlab='',ylab='Craddock')

#

#Homogénéity Analysis

names(djiclim)

library(ggpubr)

ggqqplot(djiclim$Rain_Rrr)

#

library("car")

qqPlot(djiclim$Rain_Rrr)

#

#shapiro.test(rnorm(100, mean = 5, sd = 3))

shapiro.test(runif(100, min = 2, max = 4))

#require(stats)

shapiro.test(djibclim$Rain_Rrr)

#

library(tseries)

jarque.bera.test(djibclim$Rain_Rrr)

#

#Kolmogorov-Smirnov Tests

t.test(djibclim$Rain_Rrr)

wilcox.test(djibclim$Rain_Rrr)

ks.test(djibclim$Rain_Rrr)

#

#

#   2>   >>   << Program Rclim Dex

#library("RClimDex")

#rclimdex.start()

getwd()

setwd("D:/Docs.Abdi-Basid.ADAN/14.Doc.LMME/2.Project 2/Climat/RclimDex")

source("1.0.RclimDex.R.txt")

require("climdex.pcic")

require("Kendall")

require(tcltk)

#year month day without label

#I+1  i-1

#

#

#   3>   >>   << Program ClimPact2

#Download the above zip file to your computer

#(https://github.com/ARCCSS-extremes/climpact2/archive/master.zip)

# and extract it. This will create a directory named climpact2-master.

#In Windows: open R and select "File -> Change dir..."

#and select the climpact2-master directory created in step 1.

#Then type source('climpact2.GUI.r').

#(NOTE: If nothing happens, try run the additional command startss()).

#In Linux/macOS: change to the climpact2-master directory created in step 1,

#then open R in a terminal window and type source('climpact2.GUI.r').

#The first time Climpact2 is run it will install required R packages.

#This will likely require you to select a mirror from which to download.

getwd()

setwd("D:\\Docs.Abdi-Basid.ADAN\\14.Doc.LMME\\2.Project 2\\Climat\\RClimPact\\climpact2-master")

source("climpact2.GUI.r")

#

#   4>   >>   << Program RHtest

#setwd("D:\\Docs.Abdi-Basid.ADAN\\Doc.LMME\\RHtest")

#source ("RHtest_V4.r") 

#source ("RHtest.txt") 

setwd("D:\\Docs.Abdi-Basid.ADAN\\14.Doc.LMME\\RHtest\\RHtests-master")

source("RHtestsV4_20180301.r")

StartGUI()

#   4'>   >>   << Program RHtest

library(PCICt)

## Create a climdexInput object from some data already loaded in and

## ready to go.

## Parse the dates into PCICt.

tmax.dates <- as.PCICt(do.call(paste, ec.1018935.tmax[,c("year",

"jday")]), format="%Y %j", cal="gregorian")

tmin.dates <- as.PCICt(do.call(paste, ec.1018935.tmin[,c("year",

"jday")]), format="%Y %j", cal="gregorian")

prec.dates <- as.PCICt(do.call(paste, ec.1018935.prec[,c("year",

"jday")]), format="%Y %j", cal="gregorian")

## Load the data in.

ci <- climdexInput.raw(ec.1018935.tmax$MAX_TEMP,

ec.1018935.tmin$MIN_TEMP, ec.1018935.prec$ONE_DAY_PRECIPITATION,

tmax.dates, tmin.dates, prec.dates, base.range=c(1971, 2000))

## Create an annual timeseries of the CDD index.

cdd <- climdex.cdd(ci)

#

## Create a climdexInput object from some data already loaded in and

## ready to go.

## Parse the dates into PCICt.

tmax.dates <- as.PCICt(do.call(paste, ec.1018935.tmax[,c("year",

"jday")]), format="%Y %j", cal="gregorian")

tmin.dates <- as.PCICt(do.call(paste, ec.1018935.tmin[,c("year",

"jday")]), format="%Y %j", cal="gregorian")

prec.dates <- as.PCICt(do.call(paste, ec.1018935.prec[,c("year",

"jday")]), format="%Y %j", cal="gregorian")

## Load the data in.

ci <- climdexInput.raw(ec.1018935.tmax$MAX_TEMP,

ec.1018935.tmin$MIN_TEMP, ec.1018935.prec$ONE_DAY_PRECIPITATION,

tmax.dates, tmin.dates, prec.dates, base.range=c(1971, 2000))

## Create an annual timeseries of the cold spell duration index.

csdi <- climdex.csdi(ci)

#

## Create a climdexInput object from some data already loaded in and

## ready to go.

## Parse the dates into PCICt.

tmax.dates <- as.PCICt(do.call(paste, ec.1018935.tmax[,c("year",

"jday")]), format="%Y %j", cal="gregorian")

tmin.dates <- as.PCICt(do.call(paste, ec.1018935.tmin[,c("year",

"jday")]), format="%Y %j", cal="gregorian")

prec.dates <- as.PCICt(do.call(paste, ec.1018935.prec[,c("year",

"jday")]), format="%Y %j", cal="gregorian")

## Load the data in.

ci <- climdexInput.raw(ec.1018935.tmax$MAX_TEMP,

ec.1018935.tmin$MIN_TEMP, ec.1018935.prec$ONE_DAY_PRECIPITATION,

tmax.dates, tmin.dates, prec.dates, base.range=c(1971, 2000))

## Create an annual timeseries of the CWD index.

cdd <- climdex.cdd(ci)

 

#

#   5>   >>   <<  ANCLIM

 

at=read.table(file("clipboard"),header=T, sep="\t")

write(at, file = "data.txt", sep = " ")

 

 

#   6>   >>   <<  ANNEXE

pairs(datatable, col="blue", main="Scatterplots")

Y=cbind(LRY)

X=cbind(LRV, LRC, INT)

#

hist(Y, prob=TRUE, col = "blue", border = "black")

lines(density(Y))

#

OLSreg=lm(Y~X)

summary(OLSreg)

#

Qreg25=rq(Y~X, tau=0.25)

summary(Qreg25)

#

QR=rq(Y~X, tau=seq(0.2, 0.8, by=0.1))

sumQR=summary(QR)

#

anova(Qreg25, Qreg75)

#

QR=rq(Y~X, tau=seq(0.2, 0.8, by=0.1))

sumQR=summary(QR)

#

plot(sumQR)

 

 

r=read.table(file("clipboard"),header=T, sep="\t",dec=",",row.names=1)

str(r)

attach(r)

names(r)

dim(r)

 

summary(r)

 

library(Hmisc)

describe(r,num.desc=c("mean","median","var","sd","valid.n"),horizontal=TRUE)

 

sapply(r, sd)

sapply(r,range)

 

t=na.omit(r)

round(cor(t),3)

 

#e=scale(r, center=TRUE, scale=T)

library(caret)

preproc1 <- preProcess(r, method=c("center", "scale"))

n <- predict(preproc1, r);summary(n)

 

library(corrplot)

mydata.cor = cor(t, method = c("pearson"))

corrplot(mydata.cor, method = "number", type = "lower") 

#methode = number, ellipse, square, shade, color, pie,

#type = full" "upper" "lower"

 

library(RColorBrewer)

corrplot(mydata.cor, type = "upper", order = "hclust",col = brewer.pal(n = 8, name = "PuOr"))

 

palette = colorRampPalette(c("yellow", "orange", "red")) (20)

heatmap(x = mydata.cor, col = palette, symm = TRUE)

#https://cran.r-project.org/web/packages/corrplot/vignettes/corrplot-intro.html

names(t)

#http://www.sthda.com/french/wiki/parametres-graphiques

s=read.table(file("clipboard"),header=T, sep="\t",dec=",",row.names=1)

t=na.omit(s)

attach(t)

 

reg1=lm(OB~CH,data=t)

reg2=lm(OB~PE,data=t)

reg3=lm(OB~TA,data=t)

reg4=lm(OB~AR,data=t)

reg5=lm(OB~ER,data=t)

 

limx=c(0,500);limy=c(0,900)   #windows(height = 70, width = 70)

plot(-20,-20,xaxt="n",yaxt="n",xlim=limx,ylim=limy,cex.lab=1.2,xlab="Ground Based (mm/month)",ylab="Remotely sensed (mm/month)")

axis(1,cex.axis=1.3);axis(2,cex.axis=1.3)

 

par(new=TRUE)

plot(OB~CH,xlim=limx,ylim=limy,pch=20,xaxt="n",yaxt="n",col="red",xlab="",ylab="");abline(reg1,col="red",lwd=2)

par(new=TRUE)

plot(OB~PE,xlim=limx,ylim=limy,pch=20,xaxt="n",yaxt="n",col="blue",xlab="",ylab="");abline(reg2,col="blue",lwd=2)

par(new=TRUE)

plot(OB~TA,xlim=limx,ylim=limy,pch=20,xaxt="n",yaxt="n",col="brown",xlab="",ylab="");abline(reg3,col="brown",lwd=2)

par(new=TRUE)

plot(OB~AR,xlim=limx,ylim=limy,pch=20,xaxt="n",yaxt="n",col="orange",xlab="",ylab="");abline(reg4,col="yellow",lwd=2)

par(new=TRUE)

plot(OB~ER,xlim=limx,ylim=limy,pch=20,xaxt="n",yaxt="n",col="darkgreen",xlab="",ylab="");abline(reg5,col="turquoise1",lwd=2)

par(new=TRUE)

abline(0,1.33,col="black",lwd=2,lty=2)

 

legend(30,400, legend=c("CH", "PE", "TA", "AR", "ER"),horiz=T,bty ="n",

       col=c("red","blue","brown","orange","darkgreen"), lty=c(1,1,1,1,1), lwd=c(3,3,3,3,3), cex=0.8)

 

 

 

summary(reg1)

summary(reg2)

summary(reg3)

summary(reg4)

summary(reg5)

 

 

require(hydroGOF)

names(r);attach(r)

rmsec=sqrt(sum((OB-CH)^2)/length(OB));round(rmsec,3)

rmsec=sqrt(sum((OB-PE)^2)/length(OB));round(rmsec,3)

rmsec=sqrt(sum((OB-TA)^2)/length(OB));round(rmsec,3)

rmsec=sqrt(sum((OB-AR)^2)/length(OB));round(rmsec,3)

rmsec=sqrt(sum((OB-ER)^2)/length(OB));round(rmsec,3)

#rmse(sim=Data_Mon_11.325, obs=Datat_Obser)

 

require(hydroGOF)

PBIAS = 100 * (abs(sum(OB-CH)) / sum(OB));round(PBIAS,3)

PBIAS = 100 * (abs(sum(OB-CH)) / sum(OB));round(PBIAS,3)

PBIAS = 100 * (abs(sum(OB-CH)) / sum(OB));round(PBIAS,3)

PBIAS = 100 * (abs(sum(OB-CH)) / sum(OB));round(PBIAS,3)

PBIAS = 100 * (abs(sum(OB-CH)) / sum(OB));round(PBIAS,3)

#pbias(sim=CH, obs=OB, na.rm=TRUE)

 

require(hydroGOF)

# mse(CH, OB)

# mae(CH, OB)

gof(CH, OB)

ggof(sim=CH, obs=OB, ftype="dm", FUN=mean)

 

 

 

#   STDE  /  R2  /  SLOPE  /  r  /  CD / CV RSD      

################################################

r=read.table(file("clipboard"),header=T, sep="\t",dec=",",row.names=1)

 

attach(r)

names(r)

# "OB" "CH" "PE" "TA" "AR" "ER"

round(sd(PE),3)

 

reg=lm(OB~ER,data=r);summary(reg)

 

round(cor(r),3)

 

BIAS = sum(OB-TA)/length(OB);round(BIAS,3)

 

RMSE=sqrt(sum((OB-CH)^2)/length(OB));round(RMSE,3)

RMSE=sqrt(sum((OB-PE)^2)/length(OB));round(RMSE,3)

RMSE=sqrt(sum((OB-TA)^2)/length(OB));round(RMSE,3)

RMSE=sqrt(sum((OB-AR)^2)/length(OB));round(RMSE,3)

RMSE=sqrt(sum((OB-ER)^2)/length(OB));round(RMSE,3)

 

RMSD = sqrt(sum((R$OB-R$PE)^2)/length(R$OB));round(RMSD,3)

 

names(n)

#########

 [1] "JF.OB"   "JF.CH"   "JF.PE"   "JF.TA"   "JF.AR"   "JF.ER"   "MAM.OB"

 [8] "MAM.CH"  "MAM.PE"  "MAM.TA"  "MAM.AR"  "MAM.ER"  "JJAS.OB" "JJAS.CH"

[15] "JJAS.PE" "JJAS.TA" "JJAS.AR" "JJAS.ER" "OND.OB"  "OND.CH"  "OND.PE"

[22] "OND.TA"  "OND.AR"  "OND.ER"

 

 

##########

RSD = abs(sd(CH)/mean(CH));round(RSD, digits = 3)

 

 

 

 

###########################################

# display the diagram with the better model

#oldpar<-taylor.diagram(OB,CH,pch=19)

#taylor.diagram.modified(OB,CH, text="Model 1")

s=read.table(file("clipboard"),header=T, sep="\t",dec=",",row.names=1)

attach(s)

names(s)

str(s)

cbind(sapply(s, sd, na.rm=TRUE))

 

 

CHrmse=sqrt(sum((OB-CH)^2)/length(OB));round(CHrmse,3)

PErmse=sqrt(sum((OB-PE)^2)/length(OB));round(PErmse,3)

TArmse=sqrt(sum((OB-TA)^2)/length(OB));round(TArmse,3)

ARrmse=sqrt(sum((OB-AR)^2)/length(OB));round(ARrmse,3)

ERrmse=sqrt(sum((OB-ER)^2)/length(OB));round(ERrmse,3)

rbind(CHrmse,PErmse,TArmse,ARrmse,ERrmse)

 

 

OB=AIRPORT

CH=CHIRPS  

PE=PERSIANN

TA=TAMSATV3

AR=ARCV2   

ER=ERA5

 

#

par(new=TRUE)

 

require(plotrix)

oldpar<-taylor.diagram(OB,CH,add=F, pch=19,pos.cor=TRUE,

  xlab="Standard deviation",ylab="Standard Deviation",main="",

  show.gamma=TRUE,ngamma=10,gamma.col="green",sd.arcs=1,

  ref.sd=TRUE,sd.method="sample",grad.corr.lines=c(0.1,

  0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,0.95,0.99),col="blue",

  pcex=1.5,cex.axis=1.1,normalize=F)

 

# now add the worse model

taylor.diagram(OB,PE,add=TRUE,col="black", pcex=1.5,cex.axis=1.1,normalize=F)

taylor.diagram(OB,TA,add=TRUE,col="pink", pcex=1.5,cex.axis=1.1,normalize=F)

taylor.diagram(OB,AR,add=TRUE,col="brown", pcex=1.5,cex.axis=1.1,normalize=F)

taylor.diagram(OB,ER,add=TRUE,col="green", pcex=1.5,cex.axis=1.1,normalize=F)

 

 

# get approximate legend position # add a legend

legend(1.3,1.6,legend=c("JF.OB","JF.CH","JF.PE","JF.TA","JF.AR","JF.ER"),horiz=FALSE,

    pch=c(15,19,19,19,19,19),col=c("darkgreen","red","blue","brown","orange","darkgreen"), cex=0.7)

 

 

# SAISONNALITY

#par(mfrow=c(1,2))

#"red","blue","brown","orange","darkgreen"

s=read.table(file("clipboard"),header=T, sep="\t",dec=",",row.names=1)

names(s)

attach(s)

cor(s)

str(s)

 

# show the "all correlation" display

require(plotrix)

taylor.diagram(JF.OB,  JF.CH,           col="blue",       pos.cor=F,pcex=1.5,normalize=F)

taylor.diagram(JF.OB,  JF.PE,add=TRUE,  col="black",      pcex=1.5,normalize=F)

taylor.diagram(JF.OB,  JF.TA,add=TRUE,  col="pink",       pcex=1.5,normalize=F)

taylor.diagram(JF.OB,  JF.AR,add=TRUE,  col="brown",      pcex=1.5,normalize=F)

taylor.diagram(JF.OB,  JF.ER,add=TRUE,  col="green",      pcex=1.5,normalize=F)

 

require(plotrix)

taylor.diagram(MAM.OB,  MAM.CH,           col="blue",       pos.cor=F,pcex=1.5,normalize=F)

taylor.diagram(MAM.OB,  MAM.PE,add=TRUE,  col="black",      pcex=1.5,normalize=F)

taylor.diagram(MAM.OB,  MAM.TA,add=TRUE,  col="pink",       pcex=1.5,normalize=F)

taylor.diagram(MAM.OB,  MAM.AR,add=TRUE,  col="brown",      pcex=1.5,normalize=F)

taylor.diagram(MAM.OB,  MAM.ER,add=TRUE,  col="green",      pcex=1.5,normalize=F)

 

require(plotrix)

taylor.diagram(JJAS.OB,  JJAS.CH,           col="blue",       pos.cor=F,pcex=1.5,normalize=F)

taylor.diagram(JJAS.OB,  JJAS.PE,add=TRUE,  col="black",      pcex=1.5,normalize=F)

taylor.diagram(JJAS.OB,  JJAS.TA,add=TRUE,  col="pink",       pcex=1.5,normalize=F)

taylor.diagram(JJAS.OB,  JJAS.AR,add=TRUE,  col="brown",      pcex=1.5,normalize=F)

taylor.diagram(JJAS.OB,  JJAS.ER,add=TRUE,  col="green",      pcex=1.5,normalize=F)

 

require(plotrix)

taylor.diagram(OND.OB,  OND.CH,           col="blue",       pos.cor=F,pcex=1.5,normalize=F)

taylor.diagram(OND.OB,  OND.PE,add=TRUE,  col="black",      pcex=1.5,normalize=F)

taylor.diagram(OND.OB,  OND.TA,add=TRUE,  col="pink",       pcex=1.5,normalize=F)

taylor.diagram(OND.OB,  OND.AR,add=TRUE,  col="brown",      pcex=1.5,normalize=F)

taylor.diagram(OND.OB,  OND.ER,add=TRUE,  col="green",      pcex=1.5,normalize=F)

 

 

legend(30,53

,legend=c("OND.OBSERVATON","OND.CHIRPS","OND.PESERIANNCDR","OND.TAMSAT","OND.ARC","OND.ERA"),horiz=FALSE,

 pch=c(15,19,19,19,19,19),col=c("darkgreen","red","blue","brown","orange","darkgreen"), cex=0.7)

 

CHrmse=sqrt(sum((OND.OB-OND.CH)^2)/length(OND.OB))

PErmse=sqrt(sum((OND.OB-OND.PE)^2)/length(OND.OB))

TArmse=sqrt(sum((OND.OB-OND.TA)^2)/length(OND.OB))

ARrmse=sqrt(sum((OND.OB-OND.AR)^2)/length(OND.OB))

ERrmse=sqrt(sum((OND.OB-OND.ER)^2)/length(OND.OB))

round(rbind(CHrmse,PErmse,TArmse,ARrmse,ERrmse),3)

 

 

#RMS_Diff = sum(((OND.CH-mean(OND.CH))-(OND.OB-mean(OND.OB)))^2)/length(OND.OB)

OBsd=sd(OND.OB)

CHsd=sd(OND.CH)

PEsd=sd(OND.PE)

TAsd=sd(OND.TA)

ARsd=sd(OND.AR)

ERsd=sd(OND.ER)

round(rbind(CHsd,PEsd,TAsd,ARsd,ERsd),3)

 

CHBIAS = sum(OND.OB-OND.CH)/length(OND.OB)

PEBIAS = sum(OND.OB-OND.PE)/length(OND.OB)

TABIAS = sum(OND.OB-OND.TA)/length(OND.OB)

ARBIAS = sum(OND.OB-OND.AR)/length(OND.OB)

ERBIAS = sum(OND.OB-OND.ER)/length(OND.OB)

round(rbind(CHBIAS,PEBIAS,TABIAS,ARBIAS,ERBIAS),3)

 

R=read.table(file("clipboard"),header=T, sep="\t",dec=".",row.names=1)

library(caret)

preproc1 <- preProcess(R, method=c("center"))

#"center", "scale"

norm1 <- predict(preproc1,R)

print(norm1)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

require(openair)

## in the examples below, most effort goes into making some artificial data

## the function itself can be run very simply

## Not run:

## dummy model data for 2003

dat <- selectByDate(mydata, year = 2003)

dat <- data.frame(date = mydata$date, obs = mydata$nox, mod = mydata$nox)

 

## now make mod worse by adding bias and noise according to the month

## do this for 3 different models

dat <- transform(dat, month = as.numeric(format(date, "%m")))

mod1 <- transform(dat, mod = mod + 10 * month + 10 * month * rnorm(nrow(dat)),

model = "model 1")

## lag the results for mod1 to make the correlation coefficient worse

## without affecting the sd

mod1 <- transform(mod1, mod = c(mod[5:length(mod)], mod[(length(mod) - 3) :

length(mod)]))

 

## model 2

mod2 <- transform(dat, mod = mod + 7 * month + 7 * month * rnorm(nrow(dat)),

model = "model 2")

## model 3

mod3 <- transform(dat, mod = mod + 3 * month + 3 * month * rnorm(nrow(dat)),

model = "model 3")

 

mod.dat <- rbind(mod1, mod2, mod3)

 

## basic Taylor plot

 

TaylorDiagram(mod.dat, obs = "obs", mod = "mod", group = "model")

 

## Taylor plot by season

TaylorDiagram(mod.dat, obs = "obs", mod = "mod", group = "model", type = "season")

 

## now show how to evaluate model improvement (or otherwise)

mod1a <- transform(dat, mod = mod + 2 * month + 2 * month * rnorm(nrow(dat)),

model = "model 1")

mod2a <- transform(mod2, mod = mod * 1.3)

mod3a <- transform(dat, mod = mod + 10 * month + 10 * month * rnorm(nrow(dat)),

model = "model 3")

mod.dat2 <- rbind(mod1a, mod2a, mod3a)

mod.dat$mod2 <- mod.dat2$mod

 

## now we have a data frame with 3 models, 1 set of observations

## and TWO sets of model predictions (mod and mod2)

 

## do for all models

TaylorDiagram(mod.dat, obs = "obs", mod = c("mod", "mod2"), group = "model")

 

## End(Not run)

## Not run:

## all models, by season

TaylorDiagram(mod.dat, obs = "obs", mod = c("mod", "mod2"), group = "model",

type = "season")

 

## consider two groups (model/month). In this case all months are shown by model

## but are only differentiated by model.

 

TaylorDiagram(mod.dat, obs = "obs", mod = "mod", group = c("model", "month"))

 

## End(Not run)

 

 

#                ANNEXE             ####################################################3

#Taylor, K.E. (2001) Summarizing multiple aspects of model performance in a single diagram. Journal of Geophysical Research, 106: 7183-7192.

library(datasets)

library(ncdf4)

library(plotrix)

 

taylor.diagram(r,r,add=FALSE,col="red",pch=4,pos.cor=TRUE,xlab="MERRA SD (Normalised)",ylab="RCA4 runs SD (normalised)",main="Taylor Diagram",show.gamma=TRUE,ngamma=3,sd.arcs=1,ref.sd=TRUE,grad.corr.lines=c(0.2,0.4,0.6,0.8,0.9),pcex=1,cex.axis=1,normalize=TRUE,mar=c(5,4,6,6),lwd=10,font=5,lty=3)

lpos<-1.5*sd(Data_Mon_11.275)

legend(1.5,1.5,cex=1.2,pt.cex=1.2,legend=c("volcano"),pch=4,col=c("red"))

taylor.diagram(data,data,normalize=TRUE)

 

 

# fake some reference data

 ref<-rnorm(30,sd=2)

 # add a little noise

 model1<-ref+rnorm(30)/2

 # add more noise

 model2<-ref+rnorm(30)

 # display the diagram with the better model

 oldpar<-taylor.diagram(ref,model1)

 # now add the worse model

 taylor.diagram(ref,model2,add=TRUE,col="blue")

 # get approximate legend position

 lpos<-1.5*sd(ref)

 # add a legend

 legend(lpos,lpos,legend=c("Better","Worse"),pch=19,col=c("red","blue"))

 # now restore par values

 par(oldpar)

 # show the "all correlation" display

 taylor.diagram(ref,model1,pos.cor=FALSE)

 taylor.diagram(ref,model2,add=TRUE,col="blue")

 

 

 

FIN DU PROGRAMME

Enjoy It !

The Abdi-Basid Courses Institute