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

LinkedIn group page: https://www.linkedin.com/groups/ 

🏛️ tABCi Laboratories :

🧪 Abdi-Basid ADAN LABS  | Focus: Interdisciplinary Sciences & Education.

🧠 The Deep Thinking Lab | Focus: Philosophy & Human Sciences.

🌍 The EcoClimate Hub | Focus: Climate Science & Environment.

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